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A SHORT ACCOUNT OF RELAXATION METHODS 


By L. FOX (National Physical Laboratory, Teddington, Middlesex) 


[Received 4 November 1947] 


SUMMARY 


The relaxation method is a process of steadily improved approximation for the 
solution of simultaneous equations, and any problem that can be formulated in 
terms of simultaneous equations can, theoretically, be solved by this method. After 
the first section of this account, in which the physical basis leading to the title and 
vocabulary of relaxation is discussed, the method is presented as a simple mathe- 
matical technique. The solution of ordinary simultaneous equations is illustrated 
by an example, and devices are suggested for increased convergence. Application 
to the problem of vibrating structures is next described, in which Rayleigh’s 
principle and relaxation are used to obtain the frequencies and modes of vibration. 
An example is included. The use of finite differences is then demonstrated, whereby 
differential equations are replaced by finite difference equivalents, and the solution 
of a differential equation is effected by solving a set of simultaneous equations, the 
unknowns being the values of the required function at pivotal points of a range or 
ranges of integration. Attention is focused on partial differential equations of the 
second order in two variables, and practical details are illustrated by an example 
of the solution of Laplace’s equation, functional values being specified on a closed 
boundary. The application to other partial differential equations, including the 
biharmonic equation and the equations of vibration of membranes and of flat plates, 
is briefly considered, and the account ends with a short summary of problems 
already successfully solved by relaxation methods. 


1. Introduction 


Or the various indirect methods for solving algebraic linear simultaneous 
equations, the method of relaxation is the most recent, the most interest- 
ing computationally, and probably the most powerful. Though similar 
in many respects to the iterative methods invented by Gauss and Seidel, 
and refined by Morris and others, it is very much more flexible. A measure 
of skill in securing rapid convergence can be developed by practice and 
intelligence to a much greater degree than in the standard methods, and 
this psychological aspect is very important to the computer, engaged for 
the most part in routine calculations. 

Any problem which can be reduced to the solution of a set of simul- 
taneous equations can in theory be solved by relaxation. Apart from the 
ordinary simultaneous equations which arise in everyday work, there are 
two other important problems which can be reduced immediately to this 
form. The first involves the determination of the frequencies and modes 


of vibration of engineering structures, and the second involves the solution 
5092.3 s 
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of a certain important class of differential equations. Examples of these 





“ee : ae all 
applications are given in §§3 and 4. . 
Most of the available literature is due to R. V. Southwell, who, with us 
a small and variable team of research workers, has been engaged since 
about 1936 with the development of the method and its application to 
a large number of important problems in engineering and mathematical 
physics. These researches have been collected in two books (refs. 1 | 
and 2), and a third embodying the latest and most difficult applications He 
is in course of preparation. Southwell first used the method in the problem at 
of the determination of stresses in loaded frameworks, and all his subse- at | 
quent work keeps the vocabulary and notation appropriate to this problem. sid 
Every set of equations is considered from this angle. If the equations arise on 
from a differential equation in two variables the framework becomes a dif 
tensioned net, and the ‘equilibrium of the joints of the loaded framework’ exe 
becomes the ‘equilibrium of the nodes of the tensioned net’. ap] 
In this account the method of relaxation will be presented as a simple to: 
mathematical technique, but for a proper appreciation of Southwell’s books foll 
it is desirable to have an understanding of the physical picture which led phy 
him to use the method for frameworks, and from which the name ‘relaxa- | 
tion’ is derived. so 
2. Algebraic linear simultaneous equations i 
a 
A. The framework analogy k 
3. § 


Suppose a framework, initially unloaded, has constraints of some sort 
applied to its joints, preventing any displacement.t If the external loads } , 5 
are now applied the framework will remain unstressed and all the load is t 
borne by the constraints, some particular constraint taking the largest . 


) V 
share of the load. Now let this constraint be ‘relaxed’, that is, allow the ‘i 
joint it controls to displace so that it will relieve the constraint of its load. o 
Some of the latter will be transferred to the framework, and some to other 
constraints. If at every stage that constraint taking the largest load is I 
relaxed, it is physically evident that more and more of the external load A 
will be transferred to the framework, and that ultimately the residual equ 
force still borne by the constraints will be less than the uncertainty of the 
applied loads. The process can now cease, the problem being solved with 
+ In this section ‘displacement’ is understood to include ‘rotation’, and ‘force’ includes 
‘couple’. The following account suggests that each joint has only one degree of freedom, It i 
that is, its movement is completely specified by one coordinate x. In general, of course, neo 
each joint has six degrees of freedom, three of displacement and three of rotation, and 5 
six coordinates are necessary for its complete definition. However, this restriction to one I 


coordinate in no way affects the ideas underlying this physical concept of relaxation. 
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hese all the accuracy attainable, having regard to the engineer’s uncertainty 
as to the exact magnitude of the applied loads. 
with This idea of the ‘systematic relaxation of constraints’, or ‘relaxation 
mae methods’ for short, can be carried out in computation. A typical frame- 
m to work equation has the form 
itical 
fs. 1 Ay, Xp + Ayo %y+...+Ay,X, = b+ Ry. (1) 
tions Here x, is the displacement of the rth joint, b, is the external load applied 
blem at this joint, and a,. is an ‘influence coefficient’, being the force exerted 
ubse at the 7th joint due to a unit displacement of the sth joint. The left-hand 
blem. side of equation (1) thus represents the force exerted by the framework 
arise on joint 1, due to joint displacements 2,, 2y,..., 2,. The residual R,, the 
nes a difference between the force exerted by the framework and the force 
work exerted on the framework, represents the load taken by the constraint 
applied to joint 1. The problem is to reduce the residuals systematically 
imple to negligible amounts by making appropriate changes in the 2’s. In the 
books following table the computational steps on the right correspond to the 
‘h led physical processes on the left. 
slaxa ' 
raBLE | 
Apply constraints and external loads. 1. Take all x 0. Initial residuals are 
R, b,. 
2. Take constraint supporting largest load 2. Find largest R,. Make a change of 
and relax it, relieving it completely of its R,/a,,in x, R, becomes zero. 
vad. 
' - Some of load is transferred to adjacent 3. Other residuals R, are altered by 
> o “ 
mstraints amounts a@,,2,. 
loads } Relax constraints systematically, always 4. Continue making changes in the 2’s, so 
oad is taking that supporting largest load. that at each stage the largest residual 
is brought to zero. 
argest z ; = 
When all constraints are relieved of load, 5. When all the residuals are negligibly 
ww the measure the final displacements, thus small, accumulate the contributions to 
» load obtaining the equilibrium position. the 2’s, thus obtaining the solution to 
3 load. I 
the equations. 
. other 
. | , . . . 
oad is B. Example of solution of simultaneous equations 
ul load As a simple illustration consider the following three simultaneous 
sidual equations: 
of the : 102, —22X>5 a 185 R,, 
d with “te — On, ‘ ‘ 
I JX, — 6254 2%,—93 R,, (2) 
includes py 225 5X%3— 49 — Rg. 
— It is required to find values of «,, 2, and 2, for which R,, Ry, and R, are 
course, —" « « - 
ion, and negligibly small. 
n to ont It is helpful to construct an ‘operations table’, showing at a glance the 
on. 


effect on each residual of a unit change in each variable. This is shown 
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matrix of the original equations. 





























TABLE 2 











dR, dR, dR; 
x, = 1 10 3 I 
Ot, = I —2 —6 —2 
dx, I I 2 —5 





Table 3 gives details of the relaxation process. 


TABLE 3 














Row no. Xy Le Xe Ry R, Rs 
I ° fe) ° —185 —93 —49 
“ 18 e 147 —31 Accurate solution 

| one 

3 oe -24 +43 -—3 | +17 my = 11395 
4 —4 +3 9 | +13 Xo = —21°89 
5 a 2 5 13 | 3 Xg 1°75 
6 2 I rz —_ 5 
ae: } 
7 {| +14 | —22 | +2 +I +1 ~-I 





tow 1 gives the residuals corresponding to the initial guess 
1, = tt = 13 = 0. 

R, is the biggest residual, and is almost completely liquidated by a change 
of +18 in x,. Row 2 records this change and the new residuals of which 
R, is now the largest. When this is liquidated (row 3) by means of a 
change of —24 in a, R, again becomes the largest of the three residuals, 
and is in turn disposed of by a change of —4 in 2, (row 4). At this stage 
R, becomes for the first time the largest residual, and is removed by a 
change of +2 in 2, (row 5). Finally, in row 6, a change of +2 in x, reduces 
the residuals to +1, +1, and —1. In row 7 the contributions to the z's 
are accumulated, and the residuals are checked by substitution in equa- 
tions (2). The residuals are now very small, and there is no point in further 
relaxation unless a more accurate solution is required. 

Devices for the acceleration of convergence of the relaxation process 
come by experience and practice. For example, the largest residual is not 
the all-important quantity. Quicker convergence is secured if at each 
stage that residual is liquidated which requires the largest ‘displacement? 
for its liquidation. This change is the largest of the quantities R,/4,,, and 
may not correspond to the largest residual, particularly if the magnitude 
of the diagonal coefficients of the equations varies at all considerably. 

+ The word ‘displacement’, short for ‘joint displacement’ of the framework analogy, 


often used in Southwell’s publications to mean the change in the value of a variable and 
is a very convenient shorthand. 


in Table 2. The ‘matrix’ of coefficients is clearly the transpose of the 
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f the ‘Under-relaxation’, in which only a portion of a residual is removed, 
| and ‘over-relaxation’, in which the sign of a residual is deliberately 
changed, are two useful devices. If not immediately apparent from the 
operations table, it is often easy to decide, as relaxation proceeds, whether 
to under-relax or over-relax. In the present example, for instance, suc- 
cessive changes in both 2, and x, are of opposite sign so that, if further 
relaxation were performed, it would pay to under-relax when liquidating 
| wesiduals R, and R,. Again it may happen, after some relaxation, that 
the remaining residuals are roughly proportional to their initial values, 
when a simple multiplying factor on the values so far obtained can be 
used with great effect. 
One of the most important devices is that of ‘group-relaxation’, in 
} which two or more variables are altered at the same time by different 
mounts (if the change in each variable is the same, the device is often 
called ‘block-relaxation’). The idea is to obtain combinations for which 
me residual is changed by a large amount, the others by relatively small 
mounts. Referring to the operations table, Table 2, two good group- 


relaxation operators would be as shown in Table 4. 


) 
TABLE 4 
R, R, Ry 
‘hang — —— 
2, : I 23 2 I 
whicl , : ° 17 
s ol 
iduals Evidently these operations for the liquidation of residuals R, and R, 
one | respectively, would be far more powerful than those of the original table. 
1 by a In practice subsidiary operations of this sort would be included as addi- 
educes tional rows in the operations table. 
the #8 Further examples of these computational tricks are given in §§3 and 4. 
dies | Meanwhile the following points may be noted. 
carcnes |. It is generally useless to make any residual exactly zero, by an exact 
hange f,/a,, in the corresponding variable. This residual will almost 
ee ertainly change again subsequently, and much unnecessary arithmetic 
Limo is been done. One of the advantages of the method lies in the fact that 
odes : easy’ operations can be used, and much of the work can be done mentally 
ment , Without undue strain. The calculating machine need be used only for the 
Aye “re ( uculation of initial residuals (if some guess other than vx, 0 were taken 
a to start with) and the checking of the final ones. 
bly. 2. In practical examples the coefficients are not usually simple integers, 
—. ut this entails little extra effort. Having calculated the residuals using the 


*xact coefficients, the latter may be rounded off to convenient numbers for 
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use in the operations table. When residuals are subsequently checked (using 
the exact coefficients) they will not agree exactly with those obtained from 
the relaxation process, but this merely requires a little extra relaxation. 

3. In common with all methods of iteration, mistakes are not serious. 
If some residual is incorrect, the mistake will be found in the final check 
and further relaxation will be necessary. This final calculation of residuals 
is of course very important, and it also pays to check them fairly frequently 
in the course of the work, say whenever they have been reduced to one- 
tenth of their previously calculated values. 

4. Again like other iterative methods, it is not necessary to start with a 
large number of figures to secure a result of given precision. Figures can be 
added as the work proceeds, and accuracy builds up, as it were, from below. 

C. Limitations of the method 

The process of liquidating at every stage the largest residual does not 
always converge: the residuals may increase indefinitely. Temple has 
proved (ref. 3) that there will always be convergence if the simultaneous 
equations are derived from the conditions that some positive definite 
quadratic function of the variables is a minimum. Framework equations 
are always of this kind, one characteristic being that the equations are 
symmetric, that is, a,. = a,,. Any set can be reduced to this form by 
normalization, that is, by obtaining a new set of equations from that 
typified by (1) in the form 


t n 
ao ~\ pe 
> RF=0 (s=1 ton). (: 
ox, 
it | 


This condition is not essential for convergence, either of relaxation or of 


standard iterative methods; in the example of §2B, which was highly 
convergent, the condition was not satisfied. Moreover, the relaxation 
method, by virtue of its greater flexibility, can be applied by a skilled 
computer to equations for which normal iterative methods would fail. 
The question of theoretical convergence is not indeed very important. 
Much more significant is the question of practical convergence, which may 
be so slow as to render the method impracticable. Equations of this sort 
are called ‘ill-conditioned’. 

Consider, for example, the following set of equations, devised by 


T. S. Wilson: 


5x,+ 1X 673-4 Dv, -23 i) 
1%,+10%,.+ 823+ Tx,—32 = 0 (4) 
6x, 8x,+ 10x54 9x,—33 = 0 
5x,+ Ta,+ 97,+107,—31 0 


the solution of which is x 
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These equations are in the normal form sufficient for convergence, but 
any indirect method of solution is quite impractical. Table 5 shows three 


sets of residuals corresponding to three distinct approximate solutions. 


TABLE 5 


; 1 R, Rs R, R, 
I orl I orl orl 
Ol oI o’ol ° I 
21 o’oo! oo! o°oo!l o*oo! 


The residuals given in the first row are quite small compared with the 
initial constants occurring in equations (4), whereas the corresponding 
:pproximate solution bears no resemblance to the correct one. Even the 
third solution, whose residuals are less than 1/10° of the initial constants, 
has errors of 10 per cent. 

Though relaxation is often impracticable in cases of this sort, especially 
with a large number of unknowns, it is hardly likely to lead to incorrect 
solutions. It is found in the course of the work that small residuals require 
very large displacements for their liquidation, and the whole process is 
tedious and slow, a sure sign of lack of condition of the equations. 

Relaxation, in fact, is most useful when the set of equations has large 
coefficients down the principal diagonal, or when sufficient effective group- 
relaxation operators can be worked out to produce a similar effect. Here 
igain the skill of the computer plays a big part, and the golden rule should 
ilways be borne in mind that any device, however revolutionary, can be 
ised, subject to the final check that the solution obtained satisfies the 
riginal equations. 

The basic principles and ideas of relaxation methods have now been 
lescribed, and it remains to demonstrate their application to two im- 


portant problems. 


3. Vibration of engineering structures 

A. Relaxation and Rayleigh’ s principle 

The first important application concerns the determination of the 
irequencies and modes of vibrating systems. In terms of simultaneous 
equations it is required to find the » non-trivial solutions of equations of 
the f, 

ne jorm 
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and the n values of A for which solutions exist. These values of A are the 
roots of the determinantal equation 


Qyj—A_ yg - = Ay =o (6) 
Qj  Agg—A . «gy 
Ain Asn . : Ann —~a 


In most cases suitable for relaxation the n roots are real and positive. 
The method of solution is a combination of relaxation methods and 
\ayleigh’s principle. Writing 


X, = 44, %,t4,.%+...+-a,,%,, ete. (7) 
equations (5) become 
X,—Ax, = 0, 
X,—Ar, = 0, ete. 


(38) 


Rayleigh’s principle states that the ‘best’ value of the frequency corre- 


sponding to a guess 2, 29,..., #,, at the mode of vibration, is given by 


* 
me 
2 yt 
A : ; (9) 
n 
Ya 
—_ 
I 


and that this value is stationary with respect to small changes in the 
modal, components. 

The modes and frequencies are obtained one at a time, the procedure 
for some particular mode being as follows: 

(i) Make a guess at the mode, and calculate the corresponding value 
of A using equation (9). 

(ii) Insert this value of A in equations (8), and calculate the residuals 


R= X.—dz, (10) 


q P 
X, having been already obtained in the calculation of 4d. 

(iii) Try to liquidate the residuals by relaxation of equations (5) (from 
which the operations table, with the calculated value of A inserted in the 
diagonal elements, is obtained in the usual way). It is not possible to 
obtain zero residuals unless the calculated value of A is a root of equation 
(6), but they can be reduced or ‘smoothed’ to smaller values than the 
original set. ‘ 

(iv) When further relaxation appears to effect no improvement in the 
residuals, take the new values 2,, calculate a new value of A and hence 
new residuals, make the necessary corrections in the diagonal elements 


of the operations table, and continue the relaxation process. 
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(vy) Continue until no further change in A takes place, and the residuals 


are effectively zero, to the accuracy required. 


To obtain a second mode and frequency it is necessary merely to make 
, different guess 2,,..., z, at a mode and repeat the process. The only 
danger is that in the course of relaxation this new guess degenerates into 
the old, and no new mode is obtained. This can be prevented by using 


the orthogonal property of normal modes, which states that if x,, %,..., % 


n 
and ¥;, Yo:--» ¥, are distinct solutions of equations (5), then 
Vt 
> Ln Yr 0. (11) 
1 


Thus the new guess z probably contains a component of the first mode z, 
which can be removed by taking, instead of <z, 
y = 2—an. 


Then if y is to be orthogonal to xX, 


> 2,(z,—ax,) = 0, (12) 
l 
I > (z,2%,)/ > (4). (13) 
i ‘3 


The mode y then contains no component of x (though it does contain 
components of normal modes not yet obtained). If it is used in place of z, 
there is no initial danger of ‘slipping back’ into the old mode x, and 
subsequent danger can be prevented by repeating the elimination process 
tany stage. In practice the initial elimination is often sufficient. 

\n alternative and sometimes more useful method for obtaining a 
second mode and frequency is as follows. Suppose a mode x, with com- 
ponents 2, 2 ,..., x, and frequency A,, has been obtained. A new determi- 
nantal equation replacing (6) is then formed, the element a,, in the rth 
row and sth column being replaced by as where 

}, X,X, 
- | A, Ss x 
Une of the roots of this determinantal equation, corresponding to the 
known mode 2, is zero. The remaining roots and modes are the same as 
those of the original problem, and can be found in a similar manner. The 
writer has found this method to work well in cases which were proving 
ufficult by the first method 
B. Numerical « vam ple 


As a simple example consider the following problem with three degrees 


| freedom 


5—A)x, 3%. 47, 0 ) 


(14) 
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Equations (7) in this case become 
X, = +5x,—3x,+42, 
X, = —3a,+27,—22, }. (15) 
Xs +4x,—2x,+ 10x, 


The first stage in the solution is shown in Table 6. 


TABLE 6 





' ¥ , : : 7 dk, bk, ry 
; A=- 5 
2 X 6 3 12 3 a 
- 7 = OX, I C 3 4 
2 R I 8 7 oX%— = 1 3 —3 2 
4 vy 2 I 2 I Ox I } a 5 





The initial guess (row 1) is taken, for want of a better, to be 
xy Xo Xs & 

Values of X,, X,, and X,, calculated from equations (15), are given in 
row 2. The value of A, calculated from equation (9), is found to be 5-0. 
Initial residuals (X,—Az,) are then obtained (row 3). The operations 
table, obtained from equations (14) with the calculated value of A inserted, 
is given on the right of the table. Comparison of the residuals and the 
operations table shows that the most efficacious change for the reduction 
of residuals is one of —2 in x,, resulting in new residuals as shown in row 4. 

At this stage it is decided for several reasons that further relaxation 
is likely to be valueless. First, the diagonal elements of the operations 
table depend on A, which, in view of the change in the whole character 
of the mode due to the single operation performed, must now be expected 
to change considerably from its initial value. Second, the residuals are 
now all of opposite sign to the corresponding modal components, and it 
is clear from equation (10) that this is a satisfactory state of affairs, in 
that a decrease in A will improve all the residuals simultaneously. 


The new approximate mode is 2, , & 1,2 1, and the 


2 , "2 
process is repeated as shown in Table 7. 

The new value of A is found to be 3-67, causing a substantial change in 
the diagonal elements of the operations table. The last row in the opera- 
tions table is a group-relaxation operator, constructed to give a more 
effective liquidation of R,. This operator is used twice,+ in rows 4 and 5, 
and effects a substantial reduction in all the residuals. The last operation 
(row 6) is rather cunning! Not only does it have the effect of making the 
residuals all of the same sign as the corresponding modal components, 

+ A single operation, combining rows 4 and 5, would have meant heavier arithmetic (cf. 


§ 2B, point 1). All the present computations, except the final calculation of A, were per- 
formed mentally. 
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TABLE 7 








I II 
(15 A — 3°67 
) 4 4 5 
am 5R, dR, SR; 
ae Pe) = is = ne - ae 
| dx, I 1°33 3 4 
13 0097 | | 
} OX> I 3 I 67 2 
Pe | 8x3 = 1 +4 2 6°33 
| ———- 
ox, = 2, Oxy = 1 2 5°34 2°33 
} }1 $7 sa = senate = -o 
” I 0°87 8°5144 
A : 3°718 
3°71 7 3°24 2595 


so that an increase in A will improve them all, but also the magnitude of 


each residual is now roughly proportional to the modal component, giving 


en in considerable all-round benefit (if they were exactly proportional, a suitable 
e 5-0, change in A would yield an exact solution). The last part of the table shows 
itions the improved mode (row 7), resulting in a A of 3-718, and the final residuals 
erted row 9) are very small. This value of A is in fact correct to three decimal 


d the places, and the modal components to two. 
iction It is clear from this example that the general remarks of §2B apply 
‘ow 4 equally here. The main difference is that the mode is indeterminate to 
cation the extent of a constant multiplier, so that the residuals could be made 
utions negligible by taking this multiplier sufficiently small, leading ultimately 
racter to the trival and useless solution 2, = x, = ... = x, = 0. To prevent this 
ected | it is advisable to maintain at least one of the modal components at a 
Is are reasonable size, using a multiplying factor if necessary. 
und it It is not always possible to decide which particular mode has been 
rs, it btained, though in practice the general shape of the fundamental, the 

most important mode, is often known from physical considerations. If the 
d the roots are well separated the identity 

sum of roots = sum of elements of principal diagonal 

re gives valuable information 
seul [ll-conditioning manifests itself when two or more roots lie fairly close 
e together. This does cause trouble, but in many practical problems only 
ee : few of the n roots and modes are required, and it appears generally 
: a ind fortunately) that these are the ones most easily obtainable by relaxa- 
: = tion methods. Determinants of orders up to the twelfth have been suc- 
nts 


essfully solved. including the more difficult case in which the unknown A 


appears in every element, and not just in the diagonal terms of the 


equations 
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4. Differential equations 


A. Derivation of the simultaneous equations 

The second problem, the solution of differential equations, provides 
the most interesting and the most important application of relaxation 
methods. As in the well-known step-by-step methods of integration, 
numerical values of the required function are obtained at pivotal points 
of a range or ranges of integration, the distance or interval between pivotal 
points being at the discretion of the computer. The two methods of 
solution are, however, quite distinct, and there is a most important 
difference in the types of differential equation for which the two methods 
are respectively suitable. Relaxation favours those differential equations 
for which the boundary conditions are specified on a closed boundary, 
step-by-step methods prefer the open boundary type. For example, a 
second-order differential equation in one independent variable needs two 
boundary conditions for its solution. Common boundary conditions are 

(i) the function or its first derivative is specified at the end-points of 

the range of integration, or 
(ii) both the function and its first derivative are specified at one end 
of the range. 

The first case would be solved by relaxation, the second by step-by-step 
methods. 

Similarly, for partial differential equations in two independent variables, 
relaxation favours the elliptic type, for example, Poisson’s equation 

V2f = &f/ea?+erf/oy” = g(x,y), (16) 

for which common boundary conditions are that either f or its normal 
derivative éf/év is specified at every point of a closed boundary. With 
other types of differential equations, for example, the parabolic equation 


co ' 55 
— = = (heat conduction, ete.), (17) 
C 


ne 2 
a ll (vibration of strings, etc.), (18) 
ox" ote 
the boundary conditions are seldom specified on a closed boundary, and 
relaxation methods are generally useless. There are, of course, satis- 
factory step-by-step methods, both by hand and by machines such as the 
differential analyser, for the solution of such equations. 

The simultaneous equations, whose unknowns are the values of the 
required function at pivotal points of the range or ranges of integration, 
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come from the replacement of derivatives by finite-difference expressions. 
For first and second derivatives, the equations are respectively as follows: 


df Dai Dae 
h{ I \ = 6} — d+ 05 — «+. | 


dz} 6 30 

- : ; , (19) 
(—2) a SS —_ M4 — 
h zal 0 12°°! 90° 


where h is the interval of tabulation of the function f(x), and 8% is the rth 
central difference of the function at the pivotal point s, arranged according 


to the following scheme: 


P(2 
2/ J_o Ons 
oO oO 
J - Bt. 
” ee : ete. (6d) (814-4 54), etc.) 
0 ; of = ? s 
O4 O14 
ae ° 83 8 
re) 03 
2h i> 5 ; of 


Any of the differences occurring in equations (19) can be expressed in 
terms of functional values. Treating the first and second differences only 
in this way, the equations become 

is 1 P 1 $3 1 §$5 
hf 2\ I J 1)— 6 %OT 30 09 «++ | 
ar” ; of 1541 1 96 . 
h*fo (fi +-Sa—2fo)—iz %0+ 90% --- ) 


Many important problems involve the Laplacian operator 


(20) 


> 9 


; C o- : . ‘ , 
V? -+.—. in plane Cartesian coordinates x, y. 
Cx" cy* 
For problems of rotational symmetry, using cylindrical polar coordinates 
r, 0,2), we have ° ; oe 
oo 28 o 


V2 ' 
9 | 


or? | r Or | G22 
With the aid of equations (20), and with reference to Fig. 1, these two 
operators can be expressed in finite-difference form as 


h?V2, f, (fittadet(hthady—sfot+A: | 


aes . sy = (21) 
where A 5(06 .74 86, y) +-36(8§. -4 09, y) a 
and 
R22 h . h . . . . 
“Vs Ig = f,(2 1 Dp +f 1 —~s + (fit vz—4fot+A, 
where : aw ), (22) 
h ; : 7 P ee 
A 4 55; + 36 00,, ove) (96, és 86,-) T 90(50,r+ 99,2) — ae 
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In these equations the subscripts 2, r, ete., mean ‘in the x-direction’, ‘in 
the r-direction’, etc., and A, a function of differences higher than the 
second in the two directions, will be called ‘the difference correction’. 
For convenience the interval of tabulation A is taken to be the same in 
both directions. 








lin or F 
h 
oo 
xorz 
- hk = 
=| 
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Now imagine the region over which the solution is required, and on the 
boundary of which some condition is specified, to be divided into a square 
mesh by lines parallel to the axes. Ignoring for the moment the difference- 
correction A, equations like (21) and (22) can be used to provide, at each 
mesh point, an equation connecting the values of the required function 
at this and surrounding points. If there are n mesh points there will be 
n simultaneous equations, and the solution of these by relaxation will 
provide an approximate solution to the differential equation. The number 
of equations may be large, but the number of unknowns present in each 
equation is small, a situation ideal for relaxation methods. 


B. Numerical example 
Consider, for example, the solution of Laplace’s equation 


af, ef 


¢ i oe 9 
Ox*  oy* 





= 0 (23) 


over a rectangular area, on the boundary of which values of f are given. 
The relevant finite-difference equation, obtained from equation (21), is 
(fi +f tlh +f-s)y—“fe = 0. (24) 
Now consider Fig. 2, in which equation (24) is to be satisfied at every 
internal mesh point, boundary points having an assigned value. It is clear 
from the diagram that the value at every internal point is an unknown, 
but only the points not adjacent to the boundary, that is, those marked 
with a cross, are surrounded by four unknown values. For points marked 
with an open circle, one of the terms in equation (24) is a known boundary 
value, and for the points marked with a closed circle, two of the /’s have 
known values. Thus for points marked with a cross there are five 
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unknowns in each equation, for points marked with an open circle there 


are four, and for points marked with a closed circle only three. 


The method of 
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follows familiar lines. 


Plausible values are 


attached to each mesh point, and initial residuals calculated from the 


equatic ns 


> 
R, 


Chi 


I i vet(h +S-sahy ct 4fy. 


(25) 


The operations table is simple: clearly a change of +1 in f, alters Ry by 


+, and changes the residuals at adjacent internal points in the (2, y)- 

















directions by 
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Fic. 3. 


| (since the coefficients of f, in the expressions for the 


residuals at these points are always +1). It is found convenient to record 


values of the function to the left of the nodal points, residuals to the right. 


Figs. 3 (a), (b), and (c) show the effect of a unit change in the functional 
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value at each of the three types of points of Fig. 2. The boundary is shown 
with bold lines, the mesh with fine lines. 

As an illustration of the method of setting out the work, Fig. 4 gives 
full details of the relaxation process for the solution of equations of type 
(24), for a square boundary on which functional values are given. Only 
four internal points are taken, each being of the type of Fig. 3(c). . 
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The initial value zero is attached to each internal point, and the initial 
residuals, calculated by equation (25), are recorded on the right of the 
nodal points. The relaxation process then proceeds, just as for any set 
of simultaneous equations, to impose displacements, liquidating at every 
stage the largest residual. These displacements are recorded to the left 
of the nodal points, the residuals still unliquidated to the right. At the 
stage finally reached in Fig. 4, the remaining residuals show that no further 
improvement is possible, to this order of precision. The displacements are 
then accumulated, and the final residuals checked by equation (25). The 
result is shown in Fig. 5. Evidently no errors were made in the relaxation! 

Many of the remarks made in §2B apply equally here. For example, 
only occasionally in Fig. 4 is a residual liquidated completely and exactly, 
multiples of 10 being always used until the residuals are sufficiently small 
to warrant single-figure displacements. 
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The devices of under- and over-relaxation and block and group displace- 
ments require special comment. It is not always easy, in the solution of 
a general set of simultaneous equations, to decide when to over-relax or 
how to construct effective group operations. For difference equations of 
the type (24), the problem is quite easy, owing to the simplicity of the 


operations table. Clearly over-relaxation is necessary when residuals of 
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the same sign are close together, for liquidation of, say, a positive residual, 
requires a positive displacement, which imposes a positive addition to 
residuals at surrounding points. Similarly when residuals of different sign 
are close together, it pays to under-relax. In the work shown in Fig. 4, 
two over-relaxations were performed, in each case a change of -+-20 being 
made at a point at which the residual was +60. 

It is clear from Fig. 3 that the algebraic sum of the residuals changes 
!,and only if, some point adjacent to the boundary is relaxed:+ relaxation 
it any other internal point (of type 3 (a)) has no effect on the algebraic 
sum. Consequently residuals predominantly of one sign can only be made 
‘0 vanish by ‘pushing’ them towards the boundary, and this is often done 
most effectively by block-relaxation. (It is important to ensure the 


It will be noticed that the word ‘relax’ is loosely used in many contexts. Thus we 


“aX , amongst other things, a problem, an equation, a residual, a constraint, and now 
point! 


5092.3 
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effective vanishing of the algebraic sum of residuals, since small residuals 
of one sign may in the aggregate require quite large displacements for 
their liquidation.) 

A typical block-relaxation is shown in Fig. 6. 
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Fic. 6. 


A unit positive displacement applied at every internal point (the 
boundary is shown in bold lines) produces changes in the residuals as 
shown in the figure, giving a total change of —12 in the algebraic sum. 

Now consider the problem shown in Fig. 7, where the boundary values 
are specified as shown, and the initial value zero attached to each mesh 
point results in residuals all of the same sign, with an algebraic sum of 2400. 

Application of the block-relaxation operator of Fig. 6, with 200 as unit 
instead of 1, will clearly reduce the algebraic sum to zero, and leave the 
position shown in Fig. 8. 

The residuals can now be liquidated quickly by point relaxation, and in 
fact displacements of +50 at the points with residuals +200 in Fig. 8, 
and —50 at the points with residuals —200, reduces all the residuals to 
zero and completes the solution. 

Another useful device is a group relaxation of a particular kind, of great 
value in the solution of difference equations of the type 


(fAittavdet(fitt1),—4fo = constant, K say, (26) 


where f is zero on the boundary of the region considered. At some stage 
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in the relaxation it may arise that the residuals are predominantly of one 


of 


sign, though no one is very large. The device consists in calculating the nu 


average value of the left-hand side of equation (26), and then multiplying as 


the solution so far obtained by a constant such that this average value me 


is equal to the required value K. Clearly the algebraic sum of the remain- | og) 
ing residuals will be small, positive and negative signs mingling freely, 
and point relaxation should then be effective. 

It is not possible in this account to list all the tricks which a computer 
acquires for the acceleration of convergence of the relaxation process. 
As in all other branches of numerical mathematics, facility in the relaxa- 
tion technique comes with practice, and the need for such devices becomes 
apparent only when one tries to solve a differential equation without 
using them. 


C. Special problems in the solution of second-order differential equations 
The example given in the last section related to the solution of the 
easiest type of partial differential equation. The operations table was 
simple, and was effectively the same at all points. A more general case 
of a linear second-order partial differential equation might lead to a finite- 
difference equation of the form 
(Af,;+ Bf. vat (Cf, + Df. iy—F, 0 = F, (27) 
where the coefficients A,..., F are either constants or functions of the 
independent variables x and y. There is clearly no new principle involved 


in solving by relaxation this more general problem, and no further comment 

will be made. Similarly the application of the method to differential 

equations in one independent variable is straightforward and obvious. 
There is, however, a problem in the solution of differential equations 


fica 


additional to that of relaxing simultaneous equations. The latter can be 
poi 


performed to any required degree of accuracy, but the question remains 
to what extent the simultaneous equations are an adequate representation vas 
of the differential equation. The chief sources of error lie in (i) the neglect 
of the difference correction A, and (ii) the derivation of difference equa- lort 


tions for use near a curved boundary. These problems are discussed only 

briefly here; a fuller treatment, with examples, is contained in ref. (9). | whe 
In most published work the difference correction is frankly ignored, and 

the mesh made sufficiently fine to ensure that the neglect is not serious. 

The writer prefers to use as coarse a mesh as possible, on which a first 


approximate solution is obtained by neglecting the difference correction. } Fro 
The latter is then calculated at each mesh point from the approximate _ ont 
solution, and is inserted as a new residual, to be relaxed in turn. This 1s | _ first 


repeated until the full finite-difference equation is satisfied, the number 
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of repetitions necessary rarely exceeding two. By this means the large 
number of equations required by the first method is avoided, and doubts 
as to the accuracy of the solution are to a large extent dispelled. The 
method is particularly advantageous whenever the differential equation 
contains a first derivative, for example, in equation (22). 


The problem presented by a curved boundary is‘illustrated in Fig. 9. 


O 








ae 























A finite-difference equation of type (24) cannot be used without modi- 
ieation at a point such as 1, since it involves a knowledge of the external 
point O. One method of dealing with this problem is to express the value 


it 0 in terms of the known boundary value B and the internal points 


1,2, 3,.... The expression comes from the Gregory-Newton interpolation 
formula 
, U\ a2. ([7\a3 92 
Ie Jo xy T (3) 8 T (;) OT sees (28) 
where . 
a OB/h, 
Ao = fifo: 


As = i, 2fit+fo, ete. 


From this equation f, can be obtained in terms of as many internal points 
a as desired. For simplicity Southwell retains only the 


irst lorward difference A, in equation (28), that is, he assumes the function 


0 


on the line 1. 2. 3 


be linear at the boundary. The writer preserves the accuracy obtainable 
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by the difference correction method by retaining as many differences as 
are significant in equation (28). The resulting formulae for points adjacent 
to the boundary are then somewhat complicated, and may contain more 
than the five unknowns present in the equation for an ordinary internal 
point, but again there is no inherent difficulty in the relaxation. 


the function, instead of the function itself, is specified at the boundary. 
If the boundary is curved, the setting up of accurate finite difference 
equations is particularly difficult. One method is given in Southwell’s 
book (ref. 2), a second in ref. (4). 


5. Further problems in differential equations 

A. The biharmonic and allied equations 

The application of relaxation methods to the solution of equations 
containing the biharmonic operator 


‘ o2 C2 2 
Vi=(— += (29) 
Ga? © Oy? 


is described in the paper given as ref. (5). The principle is the same, 
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example, here connects 13 points, and is given, with the notation of Fig. 
10, by 1 


4 8 
hAV4f, = 20f,—8 > f,+2 Lie > Sr (9 
1 5 


of As 


though the technique is more difficult. The finite-difference equation, 10° , 





A problem of a similar kind is involved when the normal derivative of 
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ices as | Two boundary conditions are necessary for a unique solution, usually the 
jacent function and its normal derivative being specified. The chief feature of 


1 more relaxation with this operator is the slowness of convergence, making very 
\ternal | necessary the use of group- and block-relaxations. The problem in finite- 


differences of accurately satisfying of the second boundary condition 


tive of | is also troublesome. However, several important problems in elasticity 
ndary. have been solved with all the accuracy necessary for engineering 
ference applications. 

hwell’s Occasionally (in the theory of flat elastic plates) the biharmonic equa- 


tion can be replaced by two second-order equations of the form 


(eta) + Ave 0 | 
CxX\CxN cy 


| ; (31) 
uations Cc [ou cv : 
| be )4+Av ss 
cy\éx © ey 
29 where A is constant, and wu and v take specified values on the boundary. 


The approximate finite-difference equations, with the notation of Fig. 11, 




















: re given by 
1+A)(u,+Us)-4 A (uUg+U4)—2(1 + 2A Jug t+ }(U;—Ug+07—Vg) = 0 | (32) 
A(v,+v3)+(1+A)(v.+¥,)—2(1 + 2A )vo+ f(Us—UgtUz—Ug) = 0 
j 6 2 5 
3 oO | 
7 \ 8 
Fic. ll. 
Here there are two functions, wu and v, to record at every nodal point, 
und also the residuals R,, R corresponding to the two equations. When 
tion, 1 i u-displacement is made, five residuals R, are changed, and also four 
f Fig regi Thi i 
no siduals R, (by a much smaller amount). This would seem to be compli- 


ted, but in fact the relaxation is easier and much more convergent than 
nthe biharmonic equation. Examples are given in the paper of ref. (5). 
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B. Vibration problems 
[t was seen in § 3 that the modes and frequencies of vibration of engineer- 
ing structures can be obtained by relaxation methods. Similarly the 
equation of membrane vibration, given by 


V2f+Af = 90, (33) 
and the general equation of buckling of flat plates, 
c ow Ow a] Cow Ow 
Viw+A| —|P,——S—)+—|P,——S—}]| =0 (34) 
Ox\ * ox cy} cy\ ~ cy Ox 


(in which P,, P,, and S are known functions of x and y) can be solved by 
these methods. 

The technique is identical with that of §3. A guess is made at the 
required mode of vibration or buckling, a first approximation to A obtained 
by Rayleigh’s principle, and the mode improved by relaxation of the 
governing differential equation, expressed in finite-difference form. 

For equation (33), Rayleigh’s principle gives 


[ ( fV?f dady . 
[ ( f? dxdy , aa 

the finite-difference form of the differential equation is 
(fi tfadet(fitfady—(4—M) fo +A = 0, (36) 


and the orthogonal property of normal modes is expressed by the equation 
|| f-f,dedy =0 (r #8). (37) 


The difference correction in equation (36) is the same as that of equation 
(21), and it is an interesting characteristic of problems of this sort that 
neglect of A in the relaxation process leads to an incorrect value of ), 
but a substantially accurate mode. The A can then be improved by in- 
cluding the difference correction in the calculation of V2, for insertion in 
equation (35). Examples are given in ref: (9). 

Similar treatment of equation (34) enables the mode of buckling and 
the critical applied edge force to be determined. Examples of the vibra- 
tion of membranes (which has important electro-magnetic applications) 
and of flat plates are given in refs. (6) and (7) respectively. 

C. Miscellaneous problems and conclusion 

The following is a short list of problems involving the solution of partial 
differential equations, already successfully treated by relaxation methods. 

(a) Second-order equations 


(i) Conformal transformation: V2f = 0, f or éf/év given on boundary. 
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(ii) Torsion of bars, uniform section: V?f = const., f = 0 on boundary. 


iii) Torsion of circular shafts of varying section: 


given on boundary. 
iv) Incompressible fluid flow: V?f = 0, f = constant on fixed boundary. 


(v) Flow of compressible fluid through a nozzle: 


V*(xp) bV2x Q, 
: (Ce : (=) 
x =I] +1] }- 
| sa éy) J 
i sil 
y% = 1 on the nozzle wall, 


= 0 on the centre line, 


and the form of the function 


a fis known. 


9 


(vi) Membrane vibration: V2s-+-Aud 0, % = 0 on boundary. 
(h) Fourth orde , equations 
(i) Bending of flat elastic plates: ) V‘w a. 
Slow motion of viscous fluid: | w, éw/év known on boundary. 


li) Stretching of plates: 


S(t Brave 0 


— | L AV2y = 0 


CY \Ca CY 


u and v given on boundary. 


(iii) Centrifugal stresses in rotors. 
This is a very difficult problem to relax, the differential equations being 
of the form 


with boundary conditions 
C o 


eC fi 4 
Ar*cos(r, v) (") cos(r,y)] 5 + | 
CV r 


rj or r 
er : 1 Cys 
Br*sin(r, v) -_—_—— 


fT os 





where A and B are constants, and (r,v) is the angle between the direction 
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of r and the normal v to the boundary. Details of the method, and other 
problems in cylindrical coordinates, are given in ref. (8). 

(c) Problems in which one boundary is unspecified in advance 

In these problems some extra condition is necessary, the satisfaction 
of which enables the position of the free boundary to be fixed. 
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(i) Seepage of fluid through a porous wall. 
The differential equation for the pressure is 
V%p = 0. 
The boundary conditions are (Fig. 12): 
on AB, p=y, 
on BC, = = constant, 
cy 
onCD, p=y-yYp; 
on DE, p=0; 


and on the free surface AE, two conditions have to be satisfied: 


Cp cx 
p=090, and = = —. 
CV cs 
(ii) Oil pressure in journal bearing. 
The differential equation is 
C Cp) . moO — 
- I 1-+-c cos 6)8 oF | aa pect (1+ccos 6)? = —6cesin 6, 
6 | c6 | oz? 


where B? and ¢ are constants. 
The boundary conditions specify that p is zero on a closed boundary 
in the (z,@)-plane, but the position of the top boundary, @ = 7+, 38 
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other unknown. The value of « is obtained by satisfying an extra condition 
on this line, namely, at the point O, 
Op ‘ 
ain @. 
ction o6 
Many other important problems of this type are dealt with in ref. (2). 
The list given is sufficient to show the power of relaxation methods. 





QO=T +x 
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Finally, it may be remarked that the relaxation method is now being 
} used in the important work of tabulation of mathematical functions, 
particularly in the case in which a limited degree of accuracy will suffice. 
In this connexion it may be noted that any regular function of a complex 
variable has real and imaginary parts, each of which satisfies Laplace’s 
equation, so that tabulation of double entry functions of this kind becomes 

a practicable possibility by relaxation. 

With regard to the possibility of treating a given problem by relaxation, 
certain criteria must be borne in mind. Success can almost always be 
guaranteed if (i) the differential equations are linear, (ii) the boundary 
conditions are given on a closed boundary and are also linear, (iii) the 
solution to the problem is unique and free from singularity, and (iv) the 
problem involves at the most two independent variables. On the other 
hand, this method is so flexible that occasionally non-linear problems, 

ndary either in differential equation or in boundary condition, can be solved. 
+a, 18 Further, it is often possible to recognize the type of any singularity 
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present, and to remove this first analytically, leaving a non-singular 
problem for relaxation. Finally, though the labour of relaxation in three 
dimensions is prohibitively great, the future use of the new electronic 
calculating machines in this connexion.is a distinct possibility. 
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SUMMARY 

mixed A method is proposed whereby, using Fourier transforms, some three-dimensional 
etober boundary value problems are brought within the scope of normal relaxation methods. 

The method is illustrated by the solution of Poisson’s equation inside the region 
re and yunded by a right cylinder of any cross-section when values of the unknown 
\, 239 function are given on the curved surface and when either (i) values of the function 

r (ii) values of its normal gradient are specified over the flat ends. The steady 
mem- temperature inside a cube when the surface temperature is constant on one face 
sin al nd zero on the remaining faces is worked out as a detailed example. The solution 

s also sketched in the case of very long cylinders when only one end affects the 
elasti required function. 
1945), 

| 1. Introduction 

lids of : , ° 

\s developed by Southwell and his co-workers, relaxation methods have 
lution been very successful in the solution of boundary value problems involving 

! . mm . . . 

1947), two space-variables. The essentials of the method are the replacement of 

the partial differential equations and boundary conditions by their finite 
VIID difference approximations and the approximate satisfaction of these at 


r lists 


the nodal points of a regular (two-dimensional) network covering the region 

} in question. When the equations involve three independent variables the 
finite difference approximations can still be obtained easily, but it is 
difficult to see how a practical method can be devised to satisfy these 
‘pproximations at the nodal points of a three-dimensional network. 

The method of Fourier (and other) transforms has recently been applied 

tothe solution of partial differential equations. In this method the applica- 
tion of the transform reduces the number of independent variables by one 
ind so, in the case of a three-dimensional problem, brings it within the 
scope of normal relaxational technique. 

To illustrate the proposed method, we consider the solution of Poisson’s 
equation inside the region bounded by a right cylinder of any cross-section 
when the required function is specified over the curved surface of the 

cylinder and when either the function itself or its normal gradient is given 
nthe plane ends. In §2 we consider cylinders of finite length and give 
a detailed example. In §3 an outline of the solution is given for the 
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corresponding problem for very long cylinders in which only the condition 
at one end affects the unknown function. The method could doubtless 
be extended to problems governed by other differential equations. 


2. The finite cylinder 
(a) Given the values of the required function on the plane ends 
Take the axis of the cylinder as the z-axis and, for convenience,? let 
the length of the cylinder be 7. Then we have to find V from 
V°V+f(2,y,z) = 0, (1) 
with V = g(x,y,z) on the curved surface of the cylinder, (2) 


and V=h,(x,y) whenz=0, ) 


=h,(x,y) whenz=7, | 
where V? is the three-dimensional Laplace operator (67/é2?-+-é?/éy?+-é/é2) 
and f, g, h,, hy are specified{ functions of the variables indicated. 
Let V(n) be the finite Fourier sine transform of V, viz. 
V(n) = ( Vsinnzdz (n = 1,2, 3....). (4) 
0 
Integrating by parts and using (3), we find 
T ‘ 
| (6*V /éz?)sin nz dz = n{[h,—(— 1)"hy|—n?V (n). (5) 
0 
Multiplying equations (1) and (2) by sinnz, integrating with respect to 
z between 0 and z, and using (5) we have 
V2 V(n)—n2V (n)+F(n) = 0, (6) 
with 


V(n) = G(n) onthe bounding curve of the cross-section of the cylinder, 


(7) 
where§$ » 
F(n) = n{h,—(—1)"h.|+ | fsin nz dz, (8) 
7 0 
Gn) [ gsin nz dz, (9) 


ny 


0 
and V? denotes the two-dimensional Laplace operator (é?/éx?+é?/éy’). 
Thus F and G can be found in terms of the given functions f, g, h,, hs as 
functions of x, y, and » (numerical integration being used if necessary) 
and equations (6) and (7) specify a set of two-dimensional problems which 
can be solved by the usual relaxation technique. 


For cylinders of length c we write 7z/c in place of z. 


++ —+ 


These functions may be specified either analytically or numerically. 
§ It is assumed that f and g are such that the integrals in (8) and (9) exist. 
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Once two-dimensional maps of V have been found for integral values 
of n, inversion to V is given by 

y = 

2 lar 7 3 » 

(2/n 2! (n)sin nz, (10) 
which follows as an immediate consequence of the well-known theorem 
concerning Fourier series. The practicability of the method depends on 
the number of relaxation solutions of (6) and (7) with n= 1, 2, 3.... 
required to give sufficient accuracy for calculation from (10). In the simple 
example considered below, four or five relaxation maps were sufficient for 
the purpose in view. 

EXAMPLE. Steady temperature in a cube of side 7 with one face at constant 
temperature A the other faces at zero temperature Sj 
With the previous notation, f = h, = h, = 0,g = KY whenz = z,g = 0 
when x = 0, y= 0, and y= 7. Equations (8) and (9) give F(n) = 0, 


G(n) = V, | sin nz dz 0 (neven), = 2V,/n (n odd) ona = z, and G(n) = 0 
0 
onnz7v= 0, y O,Yy TT. 


Thus we require relaxation solutions to the problems 


Vi V(2¢+1)—(2¢q 1)?V (2¢+1) 0 (q 5 = See (11) 
with V(2g+1) = 2H/(2g+1) onx=7, (12) 


OQ on2z 0, y 0, y TT, | 
since, for n even, V(n) is clearly identically zero. Relaxation maps for 
j= 0, 1, 2, 3 are shown in Figs. 1-4 and, to avoid decimals, recorded 
values are those of 1000V(2q+1)/). These were obtained using a square 
mesh of side 7/8 and the ‘difference correction’ method recently published 
by Fox.t It should be noted that as q increases the liquidation process 
of the relaxation technique gets quicker and quicker. 

The temperature V at any point (a, y, z) is then found by reading off 
values of V(2qg+1) at the appropriate (x,y), multiplying by 

(2/a)sin(2q+-1)z, 
ind summing. 

An analytical solution of this problem is available—it is a particular 
se of a slightly more general problem solved by Carslaw,§ the solution being 
y — (erjnt) SS) Sinhle sin(2p+Dy sin(2a-+1)= 


an a sinhlz 2p+1 2q+1 
Pp 0g { 


where [2 (2p + 1)? + (2q + 1). 


> 


(13) 


+ For a cube of side c, we write mz/c, my/c, wz/c for x, y, z. 
L. Fox, Proc. Roy. Soc. A, 190 (1947), 31. 


§ H.S. Carslaw, The Conduction of Heat (Macmillan, 2nd ed., 1921), p. 105. 
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At the centre of the cube (x = y = z = 3n), (13) gives V/V, = 0-1668, 
while the value obtained from the diagrams is 
V/V = (2/2)[ 0-267 sin(z7/2) + 0-005 sin(37/2)] = 0-167. 

The method of solution given above, although not essential in this 
particular example, does in fact give numerical values with little, if any, 
more labour than numerical calculation from the analytical solution. For 
cylinders of irregular cross-section or cases in which any of f, g, hy, hy are 
specified numerically, or both, an orthodox solution would probably be 
out of the question—little extra difficulty would be caused in the treat- 
ment proposed here. 

(b) Given the values of the normal gradient on the plane ends 

The problem for solution is now given by equations (1) and (2) and, in 
place of (3), the conditions 


eV ez = h(x.y) whenz=0, | 


14 
=h,(x,y) whenz=7-z. | 4) 
We now take V(n) as the finite Fourier cosine transform of V, i.e. 
4 
V(n) = | Vcosnzdz (n= 0,1,2,3....). (15) 
0 
Integrating by parts and using (14), 
[ (€°V /éz*)cos nz dz = —[h, —(—1)"h,]—n?V (n). (16) 


0 
Multiplying (1) and (2) by cos nz, integrating with respect to z from 0 to 
and using (16) we obtain equations (6) and (7) where now 


F(n) = —[h,—(—1)"h,|+ | fcos nz dz, 
D (17) 


: 
G(n) = | gcosnz dz. 
0 
Values of V(n) are obtained for integral values (including 0) of n by the 
relaxation method and inversion to V is now given by 


V V (0) /2-+-(2/z7) » V(n)cos nz. (18) 


n=1 


3. The semi-infinite cylinder 


The equations for solution are now (1) and (2) with V or eV /éz specified 


over the flat end z = 0 and V, éV /éz tending to zero as z>. The 
@ 

: = P SID » 1. 

appropriate transforms in the two cases are now V(é)= | Vd 
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according as V or éV/éz is specified over z = 0. The procedure is then 
similar to that in §2 except that all integrals are now between 0 and ~, : 
Two-dimensional equations analogous to (6) and (7) are obtained and . 
relaxation solutions found for values o: € sufficiently close together to | 
allow numerical integration to give V from the inversion theorem 
x , 
V = (2/n) [ V2 dé. (19) 
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( 
" SUMMARY 

A number of methods of solving sets of linear equations and inverting matrices 
are discussed. The theory of the rounding-off errors involved is investigated for 
some of the methods. In all cases examined, including the well-known ‘Gauss 
elimination process’, it is found that the errors are normally quite moderate: no 

exponential build-up need occur. 
Included amongst the methods considered is a generalization of Choleski’s method 
which appears to have advantages over other known methods both as regards 
accuracy and convenience. This method may also be regarded as a rearrangement 


of the elimination process. 

TuIs paper contains descriptions of a number of methods for solving sets 
of linear simultaneous equations and for inverting matrices, but its main 
concern is with the theoretical limits of accuracy that may be obtained in 


| the application of these methods, due to rounding-off errors. 

j The best known method for the solution of linear equations is Gauss’s 
elimination method. This is the method almost universally taught in 
schools. It has, unfortunately, recently come into disrepute on the ground 

’ 


that rounding off will give rise to very large errors. It has, for instance, 


» been argued by Hotelling (ref. 5) that in solving a set of n equations we 


should keep » logy 


4 extra or ‘guarding’ figures. Actually, although 
examples can be constructed where as many as nlog,,2 extra figures 
would be required, these are exceptional. In the present paper the 
magnitude of the error is described in terms of quantities not considered 
in Hotelling’s analysis; from the inequalities proved here it can imme- 
diately be seen that in all normal cases the Hotelling estimate is far too 
pessimistic 

The belief that the elimination method and other ‘direct’ methods of 
solution lead to large errors has been responsible for a recent search for 
ther methods which would be free from this weakness. These were 
mainly methods of successive approximation and considerably more 
laborious than the direct ones. There now appears to be no real advantage 
i the indirect methods, except in connexion with matrices having special 
properties, for example, where the vast majority of the coefficients are 
very small, but there is at least one large one in each row. 

The writer was prompted to carry out this research largely by the 
practical work of L. Fox in applying the elimination method (ref. 2). Fox 
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found that no exponential build-up of errors such as that envisaged by 
5 
Hotelling actually occurred. In the meantime another theoretical investi- 
gation was being carried out by J. v. Neumann, who reached conclusions 
similar to those of this paper for the case of positive definite matrices, and 
communicated them to the writer at Princeton in January 1947 before the 
proofs given here were complete. These results are now published (ref. 6), 


1. Measure of work in a process 

It is convenient to have a measure of the amount of work involved in 
a computing process, even though it be a very crude one. We may count 
up the number of times that various elementary operations are applied in 
the whole process and then give them various weights. We might, for 
instance, count the number of additions, subtractions, multiplications, 
divisions, recordings of numbers, and extractions of figures from tables. 
In the case of computing with matrices most of the work consists of 
multiplications and writing down numbers, and we shall therefore only 
attempt to count the number of multiplications and recordings. For this 
purpose a reciprocation will count as a multiplication. This is purely 
formal. A division will then count as two multiplications; this seems a 
little too much, and there may be other anomalies, but on the whole 
substantial justice should be done. 


2. Solution of equations versus inversion 

Let us suppose we are given a set of linear equations Ax = b to solve. 
Here A represents a square matrix of the nth order and x and b vectors 
of the mth order. We may either treat this problem as it stands and 
attempt to find x, or we may solve the more general problem of finding 
the inverse of the matrix A, and then allow it to operate on b giving the 
required solution of the equations as x = A-'b. If we are quite certain 
that we only require the solution to the one set of equations, the former 
approach has the advantage of involving less work (about one-third the 
number of multiplications by almost all methods). If, however, we wish 
to solve a number of sets of equations with the same matrix A it is more 
convenient to work out the inverse and apply it to each of the vectors b. 
This involves, in addition, n? multiplications and n recordings for each 
vector, compared with a total of about 4n3 multiplications in an independent 
solution. There are other advantages in having an inverse. From the 
coefficients of the inverse we can see at once how sensitive the solution 
is to small changes in the coefficients of A and of b. We have, in fact, 
f(A), % — _(A+), 2, 
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This enables us to estimate the accuracy of the solution if we can judge 
the accuracy of the data, that is, of the matrix A and the vector b, and 
also enables us to correct for any small changes which we may wish to 
make in these data. 

It seems probable that with the advent of electronic computers it will 
become standard practice to find the inverse. This time has, however, not 
yet arrived and some consideration is therefore given in this paper to 
solutions without inversion. A form of compromise involving less work 
than inversion, but including some of the advantages, is also considered. 


3. Triangular resolution of a matrix 

A number of the methods for the solution of equations and, more 
particularly, for the inversion of matrices, depend on the resolution of a 
matrix into the product of two triangular matrices. Let us describe 
a matrix which has zeros above the diagonal as ‘lower triangular’ and 
one which has zeros below as ‘upper triangular’. If in addition the 
coefficients on the diagonal are unity the expressions ‘unit upper 
triangular’ and ‘unit lower triangular’ may be used. The resolution is 
essentially unique, in fact we have the following 

THEOREM ON TRIANGULAR ReEso.uTion. If the principal minors of the 
matric A are non-singular, then there is a unique unit lower triangular 
matrix L, a unique diagonal matrix D, with non-zero diagonal elements, and 
aunque unit upper triangular matrix U such that A= LDU. Similarly 
there are unique L’, D’, U’ such that A = U'D'L’. 

The kth diagonal element of D will be denoted by d,. The 1k coefficient 
of the equation A = LDU gives us 1,,d, uv, = a, and since 1,, = u,, = 1 
this determines d, to be a,, and u,;, to be a,;/d,; these choices satisfy the 
equations in question. Suppose now that we have found values of 1;;, w,,, 
with 7 < 7%) (that is, we have found the first 7)— 1 rows of Land columns of U) 
and the first 7)—1 diagonal elements d,, so that the equations arising from 
the first i,—1 rows of the equation A = LDU are satisfied; and suppose 
further that these choices are unique and d, 4 0. It will be shown how 
the next row of L and the next column of U, and the next diagonal 
element d;, 4 0 are to be chosen so as to satisfy the equations arising from 
the next row of A = LDU, and that the choice is unique. The equations 
to be satisfied in fact state 


Ui ok a 


Ding hj Uj (KS tp), 
Jo 


Ls ny, Uy, a; ~.— 2 Li 
{x 


— ‘ 
ioto “io iok 


jUj Uj, (hk < %). 


0 


the right-hand sides of these equations are entirely in terms of quantities 
already determined. When k iy the first equation is satisfied and can 
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only be satisfied by putting d;, = right-hand side, determining d,,._ The 
equations for k > i, can then be satisfied by one and only one set of values 
of u;,,, provided d;, 4 0. The equations for k < 1%, can also be satisfied by 
one and only one set of values of /;,;,, since each d,, is different from 0. The 
new diagonal element d;, is not 0 because the igth principal minor of A is 


equal to the product of the first 7, diagonal elements d,.. 


4. The elimination method 

Suppose that we wish to solve the equations Ax = b by the elimination 
method. The procedure is as follows. We first add such multiples of the 
first equation to the others that the coefficient of x, is reduced to zero in 
all of them (excepting the first). We then add multiples of the second 
equation to the later ones until the coefficient of 2, is reduced to zero, 
After n—1 steps of this nature we shall be left with a set of equations of 
the form > v,;7; = ¢;. From the equation v,,,2, = ¢, the unknown z, 

tSJj 


nn 


can then be found immediately, and by substituting it in the equation 


=c,_, we then find 2,_,, and so on until by 


Un—1n-1En—-1t Un-1n Tn 
repeated back-substitution we have found all the coefficients of the 
(originally) unknown vector x. This description of the elimination process 
is all that is required in order to apply it. We shall find it instructive, 
however, to look at it further from a number of points of view. 

(1) The process of replacing the rows of a matrix by linear combinations 


of other rows may be regarded as left-multiplication of the matrix by 


is 


another matrix, this second matrix having coefficients which describe the 
linear combinations required. Each stage of the above-described elimina- 
tion process is of this nature, so that we first convert the equations Ax = b 
into J, Ax = J,b and record J, A and J,b. We then convert them into 
J,J, Ax = J,J,b, and so on, until we finally have J,,_;...J, Ax = J,,_,...J,b. 
In accordance with the theorem on triangular resolution we may write 
J,,-1-.JS, = L-landJ,,_,...J, A = DU. The matrix DU is upper triangular 


that is, it has no coefficients other than zeros below the diagonal. The 


n— 


matrix L-! and its inverse L are lower triangular. 
(2) The matrix L can be very easily obtained from the matrices 
n—1 . 
J,,.., J,-;. We have in fact L = 14+ 5 (1—J,). The proof of this will 


. r=] 
be left to the reader. 

(3) There is no need for us to take either the equations or the unknowns 
in the order in which they are given. In other words, if P, Q represent 
permutations we may solve instead A’x’ = b’, where A’ = PAQ, b’ = Pb, 
x = Qx’. The permutations P, Q may be chosen bit by bit as we carry 
the process through. One popular method is to let Q be the identity, that is, 
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to take the variables in the order given, and to choose P so that the 
coefficients in the matrices J, do not exceed unity in absolute magnitude. 
This is always possible, and for almost all matrices gives a unique P. 
Alternatively, this variation of the method may be described by saying 
that P is chosen so that d, shall have the largest possible value, and subject 
to this, d, to be as large as possible, and so on. This procedure is called 
‘taking the largest coefficient in the column as pivot’. The diagonal 
elements d,, d,..., d, are known as the first, second,..., last pivots. There 
seems to be a definite advantage in using the largest pivot in the column 
as it is likely to have smaller proportionate errors than other possible 
pivots, and saves us from the embarrassment of getting a pivot which is 
little different from zero. It is possible that there is also a further advan- 
tage in choosing the largest coefficient in the matrix as pivot. 

(4) The leading terms of the work involved in solving a set of n equations 
by the elimination method are as follows: 4n+-O(n?) multiplications and 
recordings of which 4n?+-O(n) recordings involve the vector b. 

(5) If, after we have solved one set of equations Ax = b, we are asked 
to solve a second set Ax’ = b’ with the same matrix A, we have only to 
operate on b with the matrices J,,...,J,,_, the values of which may be 
supposed to have been kept for reference, and then solve DUx = J,,_,...J,b. 
In other words, if the matrices J,....,J,,_, have been kept (amounting to 


a 


in(n—1) numbers) the work involved in solving a second set with the 
same A is that part of the original work which involved b, namely, 
in’+-O(n) multiplications and n recordings. 

This process may also be expressed in another form, which appears to 
be quite different, but actually is an identical calculation. As mentioned 
in (2), the triangle L in the resolution A = LDU may be obtained imme- 
diately from the matrices J,,...,J,,_,. If we put DUx’ = y’ we shall then 
have Ly’ = b’. The equations Ly’ = b’ may be solved for y’ by one back- 
substitution process and then the equation DUx’ = y’ solved by a second 
back-substitution. 

(6) As we have described it, the matrices J, A, J.J, A,..., are all written 
down in full. Actually, however, we are not really interested in all the 
coefficients of all these matrices. All we need in the end are J,,..., J,,_, 
md J,,_,..J,A. It is sufficient, therefore, to calculate all coefficients of 
J,-1-J,A, and those coefficients of J,...J,A which are required for the 
letermination of J,,,. If we write A” for J,...J, A we have 


Alp = AYP +(5,);,,AG-” (i >), 
where (J.) aes ieee 


r/ir aes A(r-1)’ 
or 
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and by addition 


Aj? oat A,;;+ 2 (Ss)ia As ”- 
If i < r we have Aj) = A¥~” and so 
n—1 
Ay et Ajj;+ 2 (sis Af, 
r—1 
Ai,+ z (Js)is aS 
s=1 
(J,)ir = = r—1 ; ° 
A, +S GAS? 


Thus we can obtain the numbers actually required (A‘, (J,),,) without 
recording intermediate quantities. This variation of the elimination 
method will be seen to be identical with the method (1) of §6 (the 
‘unsymmetrical Choleski method’). 

This form of the elimination method is to be preferred to the original 
form in every way. The recording involved in the work on the matrix is 
reduced from 4n*+-O(n?) to n?+O(n), and the rounding off is at the same 
time made correspondingly less frequent. 








(7) The elimination method may be used to invert a matrix. One | 


method is to solve a succession of sets of equations Ax”) = b”, where 
b” = 6,,.. The total work involved in the inversion is then n?+0O(n?) 
multiplications. Alternatively, we may invert the matrices L and DU 


separately by back-substitution and then multiply them together. The 


work is still n3+-O(n?) multiplications. 

(8) When the matrix A is symmetric, the matrices L and U are trans- 
poses, and it is therefore unnecessary to calculate both of them. The best 
arrangement is probably to proceed as with an unsymmetrical matrix, but 
to ignore all the coefficients below the diagonal in the matrices A”. These 
coefficients are all either zero or equal to the corresponding elements of the 
transpose. This fact enables us to find the appropriate matrices J, at each 
stage. 

(9) The elimination method can be described in another, superficially 
quite unrelated form. We may combine multiplication of rows and addition 
to other rows with multiplication of columns and adding to other columns. 
In other words, we may form a product J,,_,...J, AK,...K,,_,, and try to 
arrange that it shall be diagonal. The matrix J, is to differ from unity 
only in the rth column below the diagonal, and K, is to differ from unity 
only in the rth row above the diagonal. If we carry out the multiplications 
by Jj,...,J,,_, before the multiplications by K,,..., K,,_,, then it is clear that 


te in 


we have only the elimination method, for in either case we form J, A, | 
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n- 


J, J, A,... and the multiplications by K,,..., K,_, which come after actually 


involve no computation; they merely result in replacing certain coefficients 
in the matrix J,,_,...J, A by zeros (compare note (2)). It is not quite so 
clear in the case where the order of calculation is A, J, A, J,AK,, 
J,J, AK,,.... In this case, however, the right-multiplications do not alter 
that part of the matrix which will be required later; in fact, they again 
do nothing but replace certain coefficients by zeros. So far as the subse- 
quent work is concerned, we may consider that these right-multiplications 
were omitted, and that we formed J,,_,...J, A as in the elimination method. 

When this method is used and we choose the largest pivot in the matrix, 
it is clear that all the coefficients of J, and of K, do not exceed unity. 
This provides one proof that when the largest pivot in the matrix is chosen 
the coefficients of L, U do not exceed unity (in absolute magnitude). 


5. Jordan’s method for inversion 

In §4(1) we mentioned that the elimination process could be regarded 
as the reduction of a matrix to triangular form by left-multiplication of it 
by a sequence of matrices J,,..., J,,-,- In the Jordan method we left- 
multiply the matrix A by a similar sequence of matrices. The difference 
is that with the Jordan method we aim at reducing A to a diagonal,} or 
preferably to the unit matrix, instead of merely to a triangle. 

The process consists in forming the successive matrices J, A, J, J, A.,..., 
where J, differs from the unit matrix only in the rth column, and where 


J,..J,A differs from a diagonal matrix only in the columns after the rth. 
Let us put A” = J,...J3,A, X” = J....Jj, 


we shall then have 


Ai? = AY-YP+,);- AG? (t #7), 


ij 


(J) te — (#n) 
(so that Ay = 0 if +r), 

Ay = ,),,Ag-”, 

XX? = 3,.,, 

XP = XH + (Sir XG”. 


The particular diagonal to which A is reduced is at our disposal. 
Possible choices include the following. The diagonal may be the unit 
matrix. Or we may arrange that the diagonal elements of the J, are all 
unity and tolerate the non-unit diagonal elements in J,...J,A. A third 
alternative is to arrange that the diagonal elements in J,,...J,A shall be 
between 0-1 and 1 and that the diagonal elements in J, shall be powers of 10. 


| Hereafter ‘triangle’ and ‘diagonal’ will be written for ‘triangular matrix’ and ‘diagonal 
matrix’, 
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Jordan’s method is probably the most straightforward one for inversion, 
Although it can be used for the solution of equations, it is not very 
economical for that purpose. For hand work it has the serious disadvan- 
tage that the recording is very heavy and cannot be avoided by methods 
such as that suggested in connexion with the elimination method. It may 
be the best method for use with electronic computing machinery. 


6. Other methods involving the triangular resolution 

There are several ways of obtaining the triangular resolution. When it 
has been obtained, it can be used for the solution of sets of equations, or 
for the inversion of the matrix as has been described under the elimination 
method. Possible methods of resolution are described below. 

(1) We may use the formulae given in the proof of the theorem on 
triangular resolution. This involves 4n?+ O(n?) multiplications, n?+-O(n) 
recordings. This method is closely related to Choleski’s method for sym- 
metrical matrices ((7) below), and we may therefore describe it as the 
‘unsymmetrical Choleski method’. 

(2) We may apply the elimination method, regarded as a means of 
obtaining the triangular resolution; see notes (1), (2), (6) on the elimination 
method. 

(3) We may obtain simultaneously, and bit by bit, the four triangles 
L, L-!, U, U-! and the diagonal D. The method makes use of the following 
simple facts about triangles: 

(a) If we wish to invert a triangle, but only know the values in a sub- 
triangle, we can obtain the coefficients of the inverse in the corre- 
sponding subtriangle: for example, if we know the first 5 rows of a 
lower triangle L, then we can obtain the first 5 rows of L-!. 

(b) If we know the first 7 columns of a unit lower triangle then we know 
its first r+1 rows: likewise, if we know the first 7 rows of a unit 
upper triangle we know also its first r+-1 columns. 

Let us suppose that we have carried the process to the point of knowing 
the first r rows of L, the first r—2 of L-! and r—1 of U and U-!. We 
varry on the inversion of L to obtain the (r—1)th and rth rows of L”, 
and then multiply these rows into A to obtain the rth and (r—1)th rows 
of L-1A, i.e. of DU. From this we obtain at once the rth and (r—1)th 
rows of D, and dividing obtain the rth and (r—1)th rows of U. By (b) we 
have the rth and (r+-1)th columns of U and by (a) obtain those of U~. 
Multiplying we obtain the rth and (r+-1)th columns of AU-", i.e. of LD, 
and from this the rth and (r+1)th elements of D and columns of L. By 
(b) we have the (r+-1)th and (r+ 2)th rows of L. 

We can, of course, arrange to increase r by 1 instead of 2 at each stage. 
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This is essentially Morris’s escalator method (ref. 4), so called because 
by breaking off the work at any stage we obtain the solution for one of 
the principal minors of A; the order of the minor increases in steps. 
Morris's method differs in one small point. The diagonal elements D are 
not obtained as the diagonal of L-!A or of AU-!, but by using the identity 

a;.— > (AU-!),;d7(L-!A),,,, which follows from the (kk) coefficient 
i<k 


{the matrix equation A (AU-1)D-(L-!A). 

If Morris’s method is used for the inversion of a matrix the work 
involved consists of 3n*+-O(n?) multiplications (two triangle inversions 
each 4n?+-O(n?), two multiplications of a triangle by A, each 4n3+-O(n?), 
ind one multiplication of two triangles of opposite type, 4n*+-O(n?)), and 
3n°-+O(n?) recordings (this can be slightly reduced). It does not appear 
to be especially satisfactory in either respect. 

To relate the above account to Morris’s put 
-=@,, 2 (O-) Yi {ig ON x; (E- nn, We = (Lop... - 

(4) We may look for an upper triangular matrix M such that 

M*A*AM = 1, 
that is, so that AM is orthogonal. From the first 7 rows of M (which are also 
the first r columns of M*) we can obtain the first r rows of M* because 

fits triangular character, and hence the corresponding rows of M*A* 
nd M*A*A. The equation M*A*A.M = 1 is then applied, using the 
frst r columns in the (r+1)th row of the product. This determines the 
ratios of the coefficients of M in the (r+1)th row. The (r+ -1)th diagonal 
element of the equation then determines the multiplying factor. Having 
found M and AM we obtain the inverse as M(AM)*, or we may solve 
Ax = b by forming (AM)*b and then M(AM)*b. In the terminology of 
rthogonal vectors, as described below, the formation of (AM)*b would 
be ‘expressing b in terms of the base of orthogonal vectors’. 

This method is the orthogonalization process described in ref. (3), p. 9. 
It is closely related to the Morris method for symmetrical matrices (see (5) 
below). We may apply Morris’s method by forming A*A and then looking 
ior the upper triangular matrix M to satisfy M*A*AM = 1. This would 
nly involve A through the formation of A*A and hence of MA*A. Thus 
Morris's method applied to the normalized matrix A*A differs from the 

rthogonalization process only in that M*A*A is obtained as M*(A*A) 
instead of as (M*A*)A. 


We now come to methods for symmetrical matrices. These can all be 
uade to provide methods for unsymmetrical matrices by normalizing the 
given matrix, that is, forming AA* from A. For instance, if we wish to solve 
Ax =b, we may form A*A and A*b, and then solve A*Ax = A*b by 
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one of these methods. This normalizing technique is, however, of doubtful 
value. The formation of A*A involves $n3+O(n?) multiplications, so that 
the work involved is greater with normalization than without, in the case 
of solving equations, and is no less for the case of inversion. Moreover, 
normalizing tends to make equations more ‘ill-conditioned’ (see § 8 below), 
(5) A scheme mentioned in note (8) under the elimination method. 


(6) We may apply the method (1), but we shall only need to find L and | 


D, since U = L*. As a slight variation we may find LD. 

(7) Another variation on (6) is to find LD?. This method is due to 
Choleski (ref. 1). The matrix LD? may involve some pure imaginary 
numbers, but no strictly complex ones. 


(8) Morris’s method simplifies considerably for symmetric matrices. | 


From the first r rows of L we can obtain the first r columns of L*-!, i.e. U-, 
by inverting. Left-multiplication by A gives the first 7 columns of AU-, 
i.e. of LD, and from this we obtain the first (r+ 1) rows of L. Again Morris 
obtains D differently. 


This method is identical with a variation of the orthogonalization 


method, applicable to symmetric matrices and due to L. Fox (ref. 2). 
Fox regards two vectors b and c as ‘orthogonal’ relative to A if (c, Ab) = 0 
(scalar product). Fox finds a set of vectors V,, Vo,..., V,, Which are ortho- 
gonal in this sense.. The vectors Av, may be used as a base for other 
vectors: we have in fact 

(b, v,) 


ina (v,, Av,) ” 


The solution of equations is effected by means of the formula 


b,v 
A“b= > (D, Vr) y 
- (v,, AV,) 
It is best to obtain V,, Vg,..., V,, by orthogonalizing the unit coordinate-axis 
vectors, that is, besides the vectors being orthogonal, v, is restricted to bea 


linear combination of e,, @,,..., €,,, or in other words, to have all coefficients 


n? 
after the rth equal to 0. In this case the vectors v, are the rows of L™, 
and the orthogonality relation is L-'\(AL-!*) = D. The orthogonalization 


process by which L-! is found is identical with the inversion of AL~'*D™. 
7. Measure of the magnitude of a matrix 

There are a number of ways in which the magnitude of a matrix may 
be measured by a real number. They include: 


The norm. The norm N(A) of the matrix A is given by 


N(A) = (trace A*A)! = (2 a?,\}. 
a 
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loubtfu The maximum expansion B(A). This is given by 
So that | Ax (Ax, Ax)! 


ies B(A) —- = = max = 


reOver = ° » ee » ° 

bel The maximum coefficient M(A). This is the largest coefficient in the 
CLOW - 

matrix : 
10d. M(A) = max |a;; 
i,j 
d L an ' ; : 
Of these measures one of the first two above is probably of greatest 

theoretical significance. In this paper we deal chiefly with the maximum 
aue to 


coefficient, since it is the most easily computed. 
aginary é ae ; 
A number of inequalities relating these are listed below. 


atin M(X+Y) < M(X)+M(Y) (7.1) 
eB M(XY) < nM(X)M(Y) (7.2) 
f AU- B(X-+Y) : - )+ BY) (7.3) 
: Ment B(XY) < B(X)B(Y) (7.4) 
N(X-+Y) : Pm N(Y) (7.5) 
lizatior N(XY) < N(X)N(Y) (7.6) 
‘ref. 2) N(X) < nM(X) (7.7) 
Ab) = 0 M(X) < N(X) (7.8) 
> ortho M(X) < B(X) (7.9) 
r other B(X) < niM(X) (7.10) 
B(X) < N(X) (7.11) 
N(X) < n*B(X) (7.12) 

| 8. Ill-conditioned matrices and equations 
When we come to make estimates of errors in matrix processes we shall 
ind that the chief factor limiting the accuracy that can be obtained is 
ill-conditioning’ of the matrices involved. The expression ‘ill-conditioned’ 
ste-axis s sometimes used merely as a term of abuse applicable to matrices or 


tobea | equations, but it seems most often to carry a meaning somewhat similar 


cients to that defined below. 

of L-. Consider the equations 

re l-4a- O-9y 2-7 

lization | (8.1) 
1#f)-! 0-8a-+-1-7y —1-2) 


» and form from them another set by adding one-hundredth of the first to 


the second, to give a new equation replacing the first 
ix may 


0-786x-+- 1-709y 1-173 | (8.2) 
be 8.2 
0-800x-+- 1-700y 1-200 J 


The set of equations (8.2) is fully equivalent to (8.1), but clearly if we 
ittempt to solve (8.2) by numerical methods involving rounding-off errors 
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we are almost certain to get much less accuracy than if we worked with 
equations (8.1). We should describe the equations (8.2) as an ill-conditione/ 
set, or, at any rate, as ill-conditioned compared with (8.1). It is charae. 
teristic of ill-conditioned sets of equations that small percentage errors jn 
the coefficients given may lead to large percentage errors in the solution, 
If we are required to solve the equations Ax = b, but the coefficients used 
are those of A—S instead of those of A, S being a small matrix, then, to 
first order in S, the solution obtained will be x)+A-'Sx,, where x, is the 
correct solution. We may average the effect of this over a random popuh- 
tion of matrices S, and over the coefficients in the solution and matrix, 
and we shall find the 


R.M.S. error of coefficients of solution 
R.M.S. coefficient of solution 
Das ‘ M.S. r of coefficients of 
— | N(Ayw(A-) R.M. omer coefficic nts of A 
n R.MLS. coefficient of A 
This equation suggests that we might take either N(A)N(A-!) or 


Das 7 . : ee ee , 
— N(A)N(A-!) as a measure of the degree of ill-conditioning in a matrix. 
n : 


We will adopt the latter and call : N(A)N(A-!) the N-condition number of A. 
We will also use nM(A)M(A-!) as another measure of ill-conditioning and 
sall it the M-condition number of A. There is substantial agreement 
between the two measures, though the M-number tends to be the larger, 
especially with diagonal or nearly diagonal matrices. 

It should be noted that if all the coefficients of a matrix are multiplied 
by the same factor the condition numbers are unaltered, but that if a row 
or column is multiplied by a very large or a very small number the 
condition numbers are usually increased. For instance, the matrices 

0-8 0-6 / 0-008 0-006 
(8.3) and (8.4) 
—0-6 0-3 —0-6 0-5 

have the M-condition numbers 1-28 and 128 respectively and N-condition 
numbers 1 and 50-005. This may be considered quite a satisfactory 
example of the application of the definition. In practice one will tend to 
work with the same number of figures throughout a matrix, and the small 
values in the first row of (8.4) will prejudice the accuracy obtainable, 
because of the number of significant figures available. It is certainly true 
that a trivial modification improves the conditioning, but we should 
consider that until the possibility of this modification has been observed 
and action taken, the matrix remains ill-conditioned. 
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It is often stated that ill-conditioned matrices are ones which have 
small determinants, that is, small considering the magnitudes of the 
coeficients. This statement contains a certain amount of truth. It is 
certainly the case that bad conditioning and small determinants tend to 
go together. However, the determinant may differ very greatly from the 
above-defined condition numbers as a measure of conditioning. This may 


be illustrated by the cases of the matrices 


/] 0 0 1 O 0 ‘3 ] ] : @ l 
0 01 0}: 0 | ] ; ] 1-1 g i g ] 
\0 0 O1 0 0 0-01, l ] 1-1 E & 1-01 


ill of which have the determinant 0-01, and which have the M-condition 
numbers 30, 300, 69-3, 612, respectively, and N-condition numbers 4-77, 
47-1, 33-0, 232. 

The best conditioned matrices are the orthogonal ones, which have 
Y-condition numbers of 1. Their M@-condition numbers are mostly of the 
order of magnitude of log» (for large order n). If the coefficients of a 
matrix are chosen at random from a normal population we shall get 
Y-condition numbers of the order of n* and M-condition numbers about 
log n times greater. Thus random matrices are only slightly ill-conditioned. 

The matrices which occur in practical problems are by no means random 
inthis sense. There is a very large class of problems which naturally give 
rise to highly ill-conditioned equations. Suppose, for example, that we have 
reason to believe that some function of position in two dimensions can be 
represented by a polynomial of the fourth degree and that we wish to 
determine the coefficients. To this end we measure the values of the 
function at 25 points, and so obtain 25 linear equations for the desired 
coefficients. It may well happen that we are only able to make the 
iieasurements within a small region, and this will certainly mean that 
the equations are ill-conditioned. In such a case the equations might be 
improved by a differencing procedure, but this will not necessarily be the 
case with all problems. Preconditioning of equations in this way will 
ilways require considerable liaison between the experimenter and the 

omputer, and this will limit its applicability. 


9. The classical iterative method 

Suppose that B is an approximate inverse of A. Then we can obtain 
irom it a better inverse B, by the formula B, = 2B—BAB. If we write 
E=1—AB, E, = 1—AB,, so that E and E, give a measure of the 


> 
9 


incorrectness of the two inverses: we have E, E?, so that at each 


pplication of this process the error is essentially squared. 
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The work involved in applying this method is considerable, since jt 
involves 2n* multiplications at each stage. It may be useful in cases 
where a good approximate inverse is already available, and 1—AB has 
already been calculated, but found to be a little larger than can be 
tolerated. We may then calculate B, but carry the process no farther, 
This involves n* multiplications, but since we may write B, = B-+-BE, 
the number of figures in one of the factors (viz. in E) may be kept 
small. 

A somewhat similar type of method applies for the improvement of 
solutions of sets of equations. Suppose, for example, we have to solve the 
equations Ax = b and that we have obtained a resolution A = L.DU 
(say), somewhat inaccurately. By double back-substitution we obtain 
a solution x, of L.DUx = b, which is an inaccurate solution of Ax = b, 
We may further test this solution by forming the ‘residual’ vector 
b, = b—Ax,, and if this is too large we solve Ax = b, to obtain a cor- 
rection. In this process we do not obtain ‘quadratic convergence’ but only 
convergence in geometric progression. On the other hand, the method is 
very practical because the work involved per stage is only 2n? multi- 
plications. 


10. General remarks on error estimates. The error in a reputed 

inverse 

Error estimates can be of two kinds. We may wish to know how 
accurate a certain result is, and be willing to do some additional computa- 
tion to find out. A different kind of estimate is required if we are planning 
calculations and wish to know whether a given method will lead to accurate 
results. In the former case we do not care what quantities the error is 
expressed in terms of, provided they are reasonably easily computed. 
With these estimates we wish to be absolutely sure that the error is 
within the range stated, but at the same time not to state a range which 
is very much larger than necessary. With the second type of estimate, the 
error is preferably expressed in terms of quantities whose meaning is 
sufficiently familiar that the general run of values involved may at least 
be guessed at. We are also as much interested in the statistical behaviour 
of the errors as in the maximum possible value. 

This paper is mainly concerned with estimates of the second kind, since 
those of the first kind can be quickly dismissed. Let B be a reputed inverse 


of A. To determine its accuracy we form E = 1—AB. Then in view ol 


the inequalities (7.1), (7.2), and the equation 


A-1—B B(E- E?-+...) 
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we have 
J . eo nM (B)M(E) 
M(B—A-') < > M(BE’) < > vw M(B)M(E) = - — 
\ r=1 r=1 ( i ( Ms 1— nM(E) 
which is the required error estimate. In order to apply this inequality it 
is necessary to carry out the matrix multiplication BA, involving n° 
multiplications. However, if it is intended to apply the classical iteration 
method for improving the inverse at least once, we shall have to calculate 
E in doing so, and we shall have 1—AB, = E, = E? and therefore 


M(B.—A-~) - nM(B,)M(E,) _ n?M(B,){M(E)}* 
. 1—nM(E,) 1—n?{M(E)}? 

It should be observed that this inequality is only applicable to the 
inversion of a matrix, and not to the solution of equations. It is difficult 
to determine the accuracy of the solution of a set of equations without 
inverting the matrix. This is another reason why it is preferable to treat 
inversion rather than solution of equations as a standard process. 

When making estimates of the effects of rounding-off errors we need the 
process under examination to be rather minutely described. If, forinstance, 
, product abc is to be formed, we need to know whether it is obtained as 
ib.corasa.bc. If it is obtained as ab.c we shall need to know how many 
figures are kept in ab. This may be either a definite number of decimal or 
binary places, or a definite number of significant figures, or the number 
f figures kept may be made to depend on the results of previous calcula- 
tions. Usually, however, by a trivial modification of the quantity recorded, 
these latter cases can be reduced to one of the former. 

The variety of possible detailed calculation procedures is, of course, 
vastly greater than the list of methods which we have considered, for 
these can be subdivided into numerous alternatives which appear only 
trivially different at first sight, but which may differ very seriously from 
the point of view of error estimates. We cannot here carry out the analysis 
lor more than a very few of the procedures. These have been chosen so as 
to give bounds of error which are both reasonably small and also fairly 
imple in their analytical form. We have concentrated particularly on 
ror estimates which can be expressed in terms of the matrix A and its 

iverse. In practical work the details of the procedure must be determined 

y other considerations. With any particular procedure it will usually be 
ound possible to obtain some estimate of the type proved in this paper, 

ut usually quantities such as M(L), M(D-"), etc., will be involved. 

These can be obtained conveniently as a by-product in the calculation. 

\lternatively, one may find bounds of error by calculating 1— AB as above. 

in this case the importance of the analysis which follows is to show that 


5092.3 


ood 


xX 











302 A. M. TURING 


it is probable that the error obtained will be reasonably small if a process 
is used which is somewhat similar to one of those here considered, and that 
these methods are therefore reasonable ones to use Our main purpose in 
this paper is to establish that the exponential build-up of errors need not 
occur, and this will be proved when we have found one method of inversion 
where it is absent. 


11. Rounding -off errors in Jordan’s method 
The Jordan method was described in § 5, but we have now to specify 
the details of the rounding-off and the diagonal. We shall consider the 
case where A is reduced to a unit matrix. We assume that in the calcula- 
—-1) A(r-1) 
AG-DA 


(r—1) ? 
Av. 


tion of each quantity 


(7—1) 
Ay? — 


an error of at most e is made. How this is to be secured need not be 
specified, but it is clear that the number of figures to be retained in 
A‘ ?/AY~? will have to depend on the values of the A\’~). Likewise, we 
assume that in“the calculation of 
xo» _ PAT” 
W Ag) 

an error of at most e’ is made. It is convenient to think of these errors 
as quantities deliberately added after the accurate calculation has been 
made. If the quantities added after the calculation of A”, X” are the 
matrices S,, S'. we shall have 

J,,[...{J(J, A+S,)+S,}...]+S8, = 1, 

J,,[...{Jo(J,+S,)+S,}...]+S, = &, 

where & represents the actual matrix obtained at the end of the calculation 
as the value of A-!. 


(11.1) 


The equations (11.1) give us 


A+ > X71S, = X;}, 1.3) 
1+ >X7!S; =X 


and hence = (1+ A-! > X;1S,)A | ee Za et. (11.3) 
r . rT 


The matrix X,A is the result of the first r stages of the reduction of A 
and agrees with D in the first ry columns. This fact may be expressed in 
the equation 


(X,A—1)I, = 0, (11.4) 


7 


where I, is that matrix which agrees with the unit matrix in the first ’ 


r 


columns and with the zero matrix elsewhere. It is also clear that X, differs 
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in th 


Fron 


’ Whe 
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Process | fom the unit matrix only in the first r columns; this fact may be expressed 
nd that a4: 

in the equation 
pose i (X, 1)(1—I,) 0. (11.5) 
ed not | From (11.4) and (11.5) we now find X>?; 
versiot X~1 = AI,.+1—-I,. (11.6) 


, When we ignore the second-order terms in the rounding-off errors (11.3), 


11.6) give us 


=—A-!} A-"(> X71S,)A-1+A-! > X78, 
specily r r 
ler the > {I.+A-(1—I,)}(S, A-?—S,). (11.7) 
-alcula- f 
Let us now assume that each coefficient S, is at most « and each coefficient 
fS. at most «’. From (11.7) we can estimate the error in /-measure 
M(E—A-!) < } nf1+ M(A-)}M(S, A-!—S)) 
not he q 
ined it > n{1+M(A-!)HneM(A-!)+e’} 
vise We ‘ 
n2{1+M(A-)Me’+-neM(A-)}, (11.8) 
rin B-measure, 
B(E—A-) < ¥ BI,+A-"(1—1,){ B(S, A-!)+ B(S))} 
> errors > ‘1 B(A 1) ent B(A 1) }-¢’n?} 
iS been 
von the nif1+ B(A-)¥e B(A-1)+ ’}, (11.9) 
rin N-measure, 
V(E—A-!) < ¥ N{I.+A-"(1—I1,){N(S, A-1)+ N(S))} 
LL] a , ’ 
> fr (l—r)t*N(A-!)HneN(A-!)+ ne} 
ulation 2(n+-1)8f{14+-N(A-)Ve’+eN(A-)}. (11.10) 
| If we use the relations S,I, = S/(1—I,) = 0, which follow from the 
| restrictions on the coefficients which can suffer rounding-off errors, (11.8) 
| may be improved to 
(11.2 " 
P M(E—A4) < ne’ ™"—" MA ete + <M(A-*)). (11.11) 
(11.3) | : ° 


is result is best possible in the sense that given e, «’, 1 we can find 
m of A }% S, A so that M(S,) < «, M(S;) < e’, M(A-1) = M and the error 


.ssed in | “(&—A-!), still ignoring second-order terms, is exactly 
.. ena—1) ,. Za—] 
(11 t ile M €-re€ = «eM ° 
2 3 
firs 


We may also use (11.7) to give us an estimate of the statistical error. 
the coefficients of the matrices S, geeag S, which are not obliged to be 
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0 be 8,,..., 8, in some order, and likewise let the coefficients of S; 
which are not necessarily zero be 8x,,,..., 8p. The equation (11.7) may 
then be put in the form 


> 


(S—A-}),, = 3 Ciju Sy» 


u= 


where c;;, depends only on the coefficients of A-!. Suppose that the 


rounding-off errors s,, are independent and have standard deviation o, 


P 

and zero mean, then the mean square value of (E—A-?);; is ¥ c?,, 0%. 
u=1 

Let us put o, = 7 for u < K, o, = 7’ for uw > K and the mean square 


> 


error in Ajj! becomes 7? > c?.,+7”? > c?,,. When we substitute in the 
F441 


u= 
correct values for c;;,, we shine: 
mean square error in (A-*);, 


= 2 (A )2(A-1)%; min(K, i—1)+ 7? p' (A-1)2.,(K —i)-+ 


m, 


woe > (A-})?,, min( j, m—1) + 0 ee , 


where 7 is the standard deviation and zero the mean of each coefficient 
of S,, and 7’ is the standard deviation and zero the mean of each coefficient 
of S,. 

Also 


mean square error in (A-'),, 


ih. — = ey ee oe 
< 7? M(A-1ia 2" 19) cy a-ry2 1)(n—1 + 


3 ; 2 


») 


+ 1f2| (MCA)}A 49425 =" 


The leading term in the R.M.S. error in (A-"),; is therefore at most 
(M(A-1)}2— 
v3 


The assumptions M(S,) < «, M(S;) < ¢’ in the above analysis state in 


effect that we are working to a fixed number of decimal places both in the 
reduction of the original matrix to unity and in the building up of the 
inverse. It is not easy to obtain corresponding results for the case where 


a definite number of significant figures are kept, but we may make som 
qualitative suggestions. 


The error when working with a fixed number of decimal places arose 
almost entirely from the reduction of the original matrix, and very little 
from the building up of the inverse. This, at any rate, applies for the 


inversion of ill-conditioned matrices with coefficients of moderate size. 
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i I 
Pa seees SI x piles 
| original matrix, so that if we work to the same number of significant 


However, the coefficients of the inverse are larger than those of the 
figures in both we may expect the discrepancy to disappear. The general 


M(S,) <8M(A), M(S{) < 8’M(A-), 


M(E—A“)_ , 1 \/.. 8 
so the < n3M(A)M(A-!)p1+ — Fee | 
inal” sso \ " M(A- 1 WG) 


| idea of this may be expressed by putting 
| 
| 


> 2. o2,| There still remains the factor WA) multiplying 5’. This could be removed 
M (2 


1 square} by arranging to reduce A, not to the unit matrix, 1, but to M(A).1. This 
would be a reasonable procedure in any case, though it would be more 


eet convenient to choose the nearest power of 10 to take the place of M(A). 
We see now that it is the M-condition number nM(A)M(A-!) which 
determines the magnitude of the errors when we work to a definite 
number of figures. 

In the case of positive definite, symmetric matrices it is possible to give 

1 ]) 


| more definite estimates for the case where calculation is limited to a 
| specific number of significant figures. Results of this nature have been 
efficient} obtained by J. v. Neumann and H. H. Goldstine (ref. 6). 

efficient} It is instructive to compare the estimates of error given above with 
the errors liable to arise from the inaccuracy of the original matrix. If we 
desire the inverse of A, but the figures given to us are not those of A but 
of A—S, then if we invert perfectly correctly we shall get (A—S)-" instead 
of A-!, that is, we shall make an error of (A—S)-!—A~—}, i.e. of 

(1—A-!S)-1A-1SA-!, 
“11)] If we ignore the second-order terms this is A~'SA~!. The leading terms 


|’ | in the error in the Jordan method were A-"(> (1 —I,)S,)A-? so that we 
af 


ost might say that the greater part of the error is equal to that error which 
would have been produced by an original error in the matrix of } (1—I,)S,. 
r 
. It is possible to give error estimates also for several others amongst the 
state ) methods suggested elsewhere in this paper. This is, for instance, the case 
thin the} for the elimination method. 
p of thi The elimination method in its first phase proceeds similarly to the 
se where) Jordan process, but we only attempt to reduce A to a triangle and not 
uke some} to a diagonal: also the matrix representing the complete operation in this 
irst phase is triangular. 
eS a©rost 
ue litth 12. Errors in the Gauss elimination process 
s for tht We will consider the errors in the Gauss elimination process as con- 


ate siz ‘isting of two parts, one arising from the reduction of the matrix to the 
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triangular form, and the other from the back-substitution. Of these we 
are mainly interested in the error arising from the reduction, since this js 
the part of the process which has been most criticized. We adopt the 
description of the process given in § 4, note (1), and observe that apart 
from a slight difference in the form of the matrices J,, the reduction is 
similar to the Jordan process. As in the Jordan process, we shall assume 


that we make matrix errors S,, Sg,..., S,, in the various stages of the 
reduction of A, and vector errors S,, Sb,..., S,, in the operations on b, 


Assuming there are no back-substitution errors, and ignoring the second- 
order terms in the errors we should have: 


error inx = U-'X,, >} X;1(s,—S, U-'X, b), 
r=1 


where X, = J,...J,. Now, assuming that the process has been done with 
the largest pivot chosen from each column, we shall have .W(X>1) = 1, for 
X;1= 1+ } (1—J,) as mentioned in § 4(2). Then 

s<r 


r 


error in x,, (A-1 > X>1(s,—S, A-'b)),,| 


m 


= > (A=?) nj(X> *)sn(S-)e— > (S,)j(A 1)» db, 


jkr lp 
j>k 


n?(n+1) M(A- rye? Mn+ )), 
2° * | . 


where M(s;) < «’, M(S,) < e. 


To these errors we have to add those which arise from the back-substitu- 


M(A-1)}2M(b)e, 


tion. This consists in solving the equations DUx = L~-'b, where U is unit 
upper triangular and D diagonal. We obtain z,, first and then 2, in order 


= d>(L-"b),— > (DU),;7;. 


of decreasing r by means of the formula 2, - 
w>r 

Now if we make an error of t, in the calculation of x, from the previously 
obtained coefficients of x, then we shall have solved accurately the 
equations DUx = L-'b+ Dt, that is, we shall have introduced an error of 
U-'t, or, since A = LDU, of A-'LDt. If we arrange that M|t,) < ed;' 
the greatest error in any coefficient from this source is n2M(A-)e, and 
normally much smaller than the error arising from the first part of the 
process. Furthermore, d, will normally tend to be less than 1. 

It is interesting to note the value of the error in the last pivot, that is, 
the error in the (nn) coefficient of J,,...J, A. The matrix error in J,,...J,A 
is X,, > X;1S,, that is, since X,, L-1 = DUA-!, itis DUA~ > X;'S,,. The 

- 


‘ 
(nn) coefficient is d,, ¥ a,,,(X;1S,),,, and since M(X;1!) = 1, M(S,) <« it does 


not exceed n?d,, M(A-)e in absolute magnitude, that is, the proportionate 
error in the last pivot is at most n2M(A-)e. This cannot be very large 
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unless the matrix is ill-conditioned. With worst possible conditioning we 

find an error somewhat similar to Hotelling’s estimate. The matrix error 

in J,..J,A may be written L-! ¥ X>1S,, from which we find that the 
r 


error in the last pivot cannot exceed n*eM(L-!). But since M(L) = 1 we 
find M(L-!) < 2”-1! (and equality can be attained): that this error may 
7 


actually be as great as 2”~*e may be seen by considering the inversion of 
a matrix differing only slightly from 


l 0 0 0 
l l 0 0 
] ] 1 0 
| 
on 2 


It appears then that the error in the last pivot can only be large if L~-! is 
large, and that this can only happen with ill-conditioned equations. 
Actually even then we may consider ourselves very unlucky if L- is 
large. Normally, even with ill-conditioned equations we may expect the 
iff-diagonal coefficients of L to be distributed fairly uniformly between 

1 and 1, possibly with a tendency to be near 0. Only when there is 
1 strong tendency for negative values will we find a large L-1. 


13. Errors in the unsymmetrical Choleski method 

When obtaining the triangular resolution of a matrix by the method of 
the theorem (§ 3) it is convenient to think of the process as follows. We 
are given a matrix A and the matrices L and DU (= W, say). We form 
the product LW coefficient by coefficient. When calculating any one of 
the coefficients of LW, we always find that the data are incomplete to the 
extent of one number, and we therefore choose this number so as to give 
the required coefficient in A. The unknown quantity when forming a;; is 
ilways either /,, or w;;. Regarding the process in this way suggests the 
following rule for deciding the number of figures to be retained. We 
always retain sufficient figures to give us an error of not more than « in 
the coefficient of A under consideration. In actual hand computation this 
tule is extremely simple to apply. Suppose, for example, that ¢ is 410-7 


4 
ind that we are forming the product (LW) q, i.e. ¥ 1p; w,4. We first form 
j=1 


Y/,;w;, accumulating the products in the machine. All the relevant 
l 

quantities should be available at this stage. We then set up the multi- 

plicand w,, which should also be known and ‘turn the handle’ until the 


quantity in the product register, rounded off to seven figures, first agrees 
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with the given value of a,, (which is assumed to have zeros in the eighth 
and later figures). All the figures in the multiplier register are then 
written down as the value of [,,. 

The theory of the errors in this method is peculiarly simple. The 
triangular resolution obtained is an exact resolution of a matrix A—S. 
where M(S) < e, and the resultant error in the inverse is A~1SA-!, and 
in any coefficient at most n*{M(A-")}*«. A similar procedure is appropriate 
in the inversion of the triangles L and W. When inverting W (say) we 
can arrange, by an exactly similar computing procedure, that its product 
with its reputed inverse differs from unity by at most «’ in each coefficient, 
i.e. LK = 1—S’, where M(S’) < « and K is the reputed inverse. Note 
the order in the product which is significant. Likewise we find a reputed 
inverse V for DU such that V.DU = 1—S” and M(S’) < e’. The error 
arising from using these reputed inverses is —(1—S”)-!'VK(1—S’)+-VK, 
or neglecting second-order terms, S’A-!+A-!S’. Finally, there is a 
possible source of error due to rounding off in the actual formation of the 
product VK. If this does not exceed «” in any coefficient, the error in any 
coefficient of the reputed inverse of A is in all at most 


n*e{ M(A-1)!2+ 2ne’M(A-!)+”. 


This paper is published with the permission of the Director of the 
National Physical Laboratory. 
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THE POSITION OF THE SHOCK-WAVE IN CERTAIN 
AERODYNAMIC PROBLEMS 


By M. J. LIGHTHILL (Department of Mathematics, 
The University, Manchester) 
[Received 25 November 1947] 
SUMMARY 

For the motions caused by the uniform expansion of a cylinder or sphere into still 
ir, and for steady symmetrical supersonic flow past a cone, first approximations 
we obtained to the strength and position of the shock-wave for small disturbances. 
It is also shown what are the correct approximations to the pressure on the surface. 
The work is mathematically rigorous, and provides some justification of approximate 
methods that have been used, for example, in the general theory of supersonic 


flow past bodies of revolution. 


1. Introduction 
TxE problems of compressible fluid flow which can be reduced to problems 
in one independent variable are the following: 
(i) The motion set up by a plane wall moving uniformly into still air, 
and 


(ii) That when it moves uniformly away from still air. 


(iii) Steady supersonic flow round a concave corner. 
(iv) Steady supersonic flow round a convex corner. 


\V 


Symmetrical supersonic flow past a cone. 

(vi) The motion set up by a cylinder expanding uniformly into still air. 

(vii) The motion set up by a sphere expanding uniformly into still air. 

Problems (i) and (ii) were solved by Riemann, and problems (iii) and (iv) 
by Meyer (1908). Problems (ii) and (iv) involve no shock-wave and will 
not be considered further here. The rest involve a shock-wave; and in 
problems (i) and (iii) it was found that the strength of the shock-wave 
was such as to produce, behind it, the uniform flow compatible with the 
boundary conditions. 

The following notations will be used throughout: y the adiabatic index, 
¢ the velocity potential, g the fluid speed, a the local velocity of sound, 
a) the velocity of sound in undisturbed air; and in problems (i), (vi), and 
vil) the notations: a) the velocity of the disturbing surface, a, M the 
velocity of the shock-wave; and in problems (iii) and (v) the notations: 
U the velocity of the main stream, a, the velocity of sound therein, 
M = U/a,, 8 = (M?—1)-, tan—a the angle of the corner or semi-angle 
of the cone, tan-!y the angle made by the shock-wave with the main 
stream. With these notations the results of Riemann and Meyer on 
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problems (i) and (iii) reduce, respectively, to relations between M and g, 
and between 7, M, and a. For small « these become asymptotically 


in problem (i) M—1 ~ }H(y+1)a, (1) 
and in problem (iii) 7n—B ~ Hy+1)M4(M2—1)-2a. (2) 


The analogous relations in problems (v), (vi), and (vii) are discovered below 
and set out in equations (43), (16), and (48) respectively. 

Problem (v) was discussed by Taylor and Maccoll (2), who obtained 
by numerical integration, for a number of values of the parameters, the 
velocity field and the position of the shock-wave. The same process was 
carried out by Taylor (3) for problem (vii). In this paper, analytical 
methods are applied to problems (v), (vi), and (vii) in the cases when 
the parameter « is small. 

When this is so the conclusions of the numerical work are as follows. 
In problem (v) the strength of the shock-wave tends to zero as « —> 0, but 
not with the extraordinary speed shown in problem (vii), for which it is 
of the order 10~1* for « = 0-2. The velocity field is in both problems 
adequately approximated by the linear theory, but the pressures are only 
approximated correctly if (in deducing them from the velocities) the full 
Bernoulli equation 2 2 

Ch ‘ a ao 9 

oT ee oT = (9) 

ot a y—l1 y—1 
is used, since in problem (vii) g? has the same order as ¢¢/ét on the sphere, 
and since in problem (v) the component (é¢/ér)? of q? is of the same order 
on the cone as the difference of the component (é¢/éz)? from its main- 
stream value (see Lighthill (1) for a fuller discussion of this point). If this 
is done, the discrepancies in the pressure noticed by Taylor and Maccoll (2) 
and Taylor (3) disappear. 

However, the linearized theory gives no approximation to the strength 
of the shock-wave; in fact it places the first disturbance to the flow on the 
‘Mach cone’ in problem (v) and on a cylinder, sphere expanding with the 
speed of sound in problems (vi), (vii). The velocity is, moreover, every- 
where continuous on this theory. This trouble cannot be avoided by the 
ordinary method of obtaining further approximations from a linear theory 

expanding the potential in powers of some parameter (say a) and 
equating coefficients—for it is found that the flow deduced by this method 
is still entirely within the Mach cone (etc.), and indeed the expansion does 
not converge near the Mach cone. (Of course this does not mean that the 
results given by such a method near the body may not be right.) 

There is nothing unusual in the divergence near a characteristic of an 
approximation to a non-linear differential equation. Examples are known 








(whe 


matl 
appl 
of th 

TI 
same 
the | 
base 
first. 
inte: 
fully 


~ 

h 
dist 
was 
only 
may 


The 
stre 


hole 


wh 


the 


he 
(iii 


th 


| beloy 


tained 
Ts, the 
‘SS Was 
lytical 


} when 


ollows 
0, but 
sh it is 
oblems 
re only 


he full 


sphere, 
> order 
malin 
If this 
coll 2 


rengt! 
on the 
ith the 
every 
by the 
theory 
x) and 
nethod 


mn dor Ss 


know! 























SHOCK-WAVE IN CERTAIN AERODYNAMIC PROBLEMS 311 


where the differential equation is partial) in many parts of applied 
mathematics, and it is hoped that the method of the present paper (here 
applied to ordinary differential equations) may later be extended to some 
of these examples. 

The method consists in replacing the equation by one which has the 
same linearization, is a better approximation than the linearization near 
the characteristic, and is soluble. A rigorous treatment is then possible, 
based on ideas from this rough method. The simplest case, treated here 
first, is problem (vi): in problems (v) and (vii), which are physically more 
interesting, the mathematics will not be rewritten but the results will be 


fully discussed. 


2. The uniform cylindrical wave 

In the uniform expansion of a cylinder from zero radius, we denote 
listance from the centre of the cylinder by r, and time since the cylinder 
was of zero radius by ¢t. Then all velocities and pressures must depend 
only on r/t (since there is no fundamental length or time involved) and we 


may ‘ite o,f 
may write db a-tf(x), x r (a,t). (4) 


The region of flow is a n M. Since the shock-wave is of uniform 
strength, the flow behind it is adiabatic and Bernoulli’s equation (3) 


holds, i.e. oI 


a’ ao| 1—(y L)(f af’ - + f’?)], (5) 
vhere accents denote differentiation with respect to 2x. 
The equat ion ot motion 
24 Ae ad\2 2 
wed — o¢ 4 Phas od +(2) ¢ $ (6) 


ot* cr oret or] ore 


therefore becomes 


[1—(y—1(f—af’ +4 ff" +2f’) = (a@—f' yf”. (7) 
The boundary conditions are: (i) at the surface of the cylinder 


Cher Ay X, 


hence f '(«) v; (il) at the shock-wave ¢ is continuous, hence f(.V/) 0: 
li) the relation between the velocity of the shock-wave and the fluid 
speed behind it (deduced from the Rankine-Hugoniot equations) gives 
that f’(.M) = 2(M—M-) (y+1) = e, say. Since there are three boundary 
conditions for a second order equation, they will enable (theoretically) a 
lunctional relation to be deduced between M and a, to which a first 
ipproximation will be obtained below. 


The linearized form of equation (7) is 
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The better approximation here adopted (as a basis for a rough argument) is 


[(it+(y—Df' lf’ +a-f’ = (2?@—2f')f’, (9) 
which has the same linearization and in which the coefficient of f” (which 
vanishes to the first order as x > 1) is given a better approximation there 
(where f is very small and f’? is small compared with f’). 

Equation (9) is easily solved. In fact, with f’ = y and 2-2 = z, it 


becomes ; 
ljdy ly 
|+ +l)y—-— — = --, 10 
yr Dy zjdz 22 20) 
whence dz/dy—2(y+1+y—)z = —2y-!, (11) 
and zyte-My+Dy — [ 2t-3e—2y+t dt+- constant. (12) 


uv 
Now y is small everywhere, so an adequate approximation to the integral 
in (12) is given by the first two terms of its asymptotic form 
¢ y+ DU y—2 2 (yy + ] je-*y+ Duy-1, 


Hence, approximately, for some K, 


z= 1—2y+1)y+ Ky’. (13) 
The boundary condition (i) now gives that z = a-? when y = «, and hence 
- l 
K~— as a->0O. (14) 
4 
xX 
The boundary condition (iii) gives that z = M-? when 


y =e = 2(.M—M-")/(y+1), 
and hence 
M-? = 1—4(M—M-)+ K.4(M—M-)?/(y+1)2, (15) 
which implies that 
3(y+1)? 
M—1~ SK 


as a> 0. (This result, for the uniform cylindrical wave, should be com- 


~ My t+ 1a! (16) 


4 


pared with the result, for the uniform plane wave, given in equation (1).) 
The rough arguments used above suggest a rigorous treatment as follows. 

While it is possible that the solution of equation (7) for a given « is not 

unique, it seems certain physically that a solution exists for every «, with 

y > 0 uniformly as «+0; this only is assumed here. Then M —> 1 as 

x > 0. Equation (7) can be rewritten as 

dy | 


[1—(y—D f+ (y+ lazy— yt Dy?—2?]= + 
da 


+2[1—(y—1)f + (y—1)axy—Hy—Dy’] = 0. (17) 
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Hence dy/dx is negative provided that y > 0 and 


(y+l)ay > 2?—14+(y—1)f+Hyt+])y’. (18) 
If « is small enough, this is true for x = M, y = e, f = 0 (since M > 1 as 
,> 0), and (18) holds when x < 1. Hence, if dy/dz is ever non-negative, 


there is a largest value of x for which it is so, and for this value 1 <a < M, 
y >e, and f = | y dx is negative. It is easily deduced that (18) holds for 
M 
this value, a contradiction. 
Hence y is strictly increasing (as « decreases from M to a), f is every- 
where O(y), and everything can be written as a function of y. Also 


ye die x?—(y IF (y+ Vey—4yt+ Dy? _ gay (19) 
“dy 1—(y—1) f+(y—l)ay—Hy—l)y? 


uniformly, and differentiation with respect to y gives that y*d*x/dy? = O(1) 
uniformly ; (17) is now rewritten (in imitation of (11)) as 


dz Z 2 
+2| : +00) = (- +a) (20) 
dy y y 
where g(y) = (2y—4a)[1—(y—1) f+(y—lD)zy—H(y—1)y?]-3, (21) 
and h(y) [1—(y—1) f+ (y+ lay—Hy t+ 1)y?|-3. (22) 
y 
It follows from (20) that, if g,(y) = exp | g(y) dy, a function which is 


evidently 1+-O(y) uniformly, then 
u 
9 


zy~*g,(y) = | - a< ay))y 29, (y)k(y) dy+- constant 


x 
UY 


y 9, (y)h(y) ( y~*g,(y)h'(y) dy-+ constant 
ygi(y)h(y)+y 1 gy(y)h'(y)— 


| yga(y)[h"(y) +9(y)h'(y)] dy +-K. (23) 


Now 

h'(y) = [(y+l)y—(y+ l)a—2y da/dyl1—(y—D) f+ (y+ Day—Fyt+l)y?|-? 
(24) 

which is O(1) uniformly; and 

h"(y) = 2[(y+1)y—(y+1l)a—2y dx/dy|? > 
[1—(y—1) f+ (y+ Day—Hyt Dy*} 3+ 
[(y+1)—(y+3) da/dy—2y d?x/dy*|> 
< (I= (y— ft (yt Dey Hy t+ Dy? }* 
(y+3) dx/dy—2y d*x/dy?+O(1) = O(y-). (25) 
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Hence, by (23), 

z= 14O0(y)+Ky{1+O0(y)] = (1+ Ky*?)[1+ O(y)]. (26) 
The boundary condition z = a-*, y = a now gives K ~ a~‘ as before. To 
apply the condition z= M-*, y =e, with greater accuracy than the 
simple deduction from (26), « = O(a?), we must estimate the integral in 
(23) when y = ec. By (25) it is 


€ 


( [(y+3)y-? dx/dy+ 2d?x/dy?| dy+ O(log «-) 
(y+3) [ (v—1)y-* dy+ O(a-!)+ O(log e-), (27) 


integrating the first portion by parts, and for the second, using the fact 
that dx/dy = O(1) at both limits of integration. Now when y = a‘, 
Z 1+-a?+ O(a) by (26), so a 1—}a?+ O(a). Hence 
4 1 4 
| (w—1)y-* dy = | (w—1)y* dy + | (v—1)y? dy 
x x as 
= O(a-%)+ O(e—a?), (28) 
and the integral of (23), with y €, Is O(a-*)+ O(e-!a?)+ O(log e-). 
Hence (23) becomes, for this value of y, 
M~* = h(e)+h'(e)+ O| (M—1)?a-4]+-O[ (M—1)a?]+ 
i 6 
+ O| (M—1)*log( M—1)-*]-+ a . pe! —1)7[1+0(M—1)], (29) 
dae i ; 
i.e. (using (22) and (24) to estimate h(e) and h'(e)) 
6(M—1){1+ O(a?) + O[(.M—1)log( M—1)-*}} 
16 
aA4 
(y+1)° 


which gives again (16), this time proved rigorously. 


= q-~4 


VM —1)*{1+O(a)+O(M—1)], (30) 


The method gives also the value of f on the cylinder. This is 
x x x 
y dx [wy]5,— | x dy y+ O(a4)—[1 + O(x)| 


A d J vy+yro *) 
M € € 


dy 


) 


x?(1 +-log 4a)+ O(alog a), (51) 


which agrees with the value given by linearized theory. Hence the 
pressure on the cylinder, if deduced from the linearized theory by equation 
(5) without the term 2 


3 


f’* omitted (since it is of the same order as af’ on the 
cylinder), will be a good approximation to the truth for small a. (The 
corresponding result for the sphere is borne out by the figures given by 


Taylor (3).) 
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3, Steady symmetrical supersonic flow past a cone 
If (r, 0, z) are cylindrical coordinates with the vertex and axis of the 
one as origin and axis, then 
db Aozj(Z), = r/2, (32) 
ince no fundamental length appears in the problem. Bernoulli’s equation 
a® = atl 1—4}(y—1){(f—af’)?+f'7}], (33) 


ind the equation of motion 


) becomes 


Vd = ayes oF od ed 4 (#)'< (34) 
Ox} cx cx cr Cxor cor or* 
comes 
(I—4(y—1){(f—af’P+f' 3 f"+atf’ +27") = {af—(1+2*)f'}f", 
(35) 


where accents denote differentiation with respect to x. The boundary 
ndition at the surface, ¢4./¢ x when r/z x, can be written 
x f(x) (1+a?)f"(«). 
At the shock-wave (2 = ») ¢ is continuous, so f(y) = U/ay. There is also 
rather complicated relation connecting f’(y) (the value of aj ' é¢/ér at 


the shhock-wave) with 7 and M, which can be approximated as 


f(y) = A(n—B)+O(n—B), (36) 

4 U/(M?—1\? 
A 8 M?—1)-! 37 
: tat we _— ' ee) 
Here » > 8 corresponds to « > 0; f may be assumed to differ only slightly 


tom Ua, for small «. The approximation obtained from (35) as (9) was 


irom (7) is therefore 
] 4 | [72 az \(f y-lf i * | "a t (y LU 4,)B(1 +.B?) f’f” 
[ U2x?/a3—2B(U /ay)(1+B2)f'|f". (38) 


hich ean be written 


| rv Udy, Q2\ ¢/ mm L os 2 
1 yt opast pays’ sty’ =o. (39) 
b= a: , Py 
1 
‘he transformation f?2~-? — z, Ua, B(1+f*)f’/aj = y, now reduces equa- 
lon (59) to equation (10): and we deduce as before equation (13). The 
undary condition at the cone is approximately f’(a) = Ua/ag, or 
M?B(1+-B?)o M*63x~ when z = B?a-*. Hence 


2 me KK. MepSq?, kK ~ M-*B-4a-4, (40) 
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as a0. But when z = f?n-?, by (36), 


y = Vay B(1+f*)A(n—B)/at+ O(n—B)? = 4(n—B)/(y+1)B+ O(n), 
(41) 
whence, by (13), 
B?n-? = 1—8(y —f)p-1+ M-§B-*&x 4 16(y— B)*(y+1)- (42) 
Hence finally, since B?7-?—1 ~ —2(y—8)p-, 
n—B ~ 3(y+1)?M8(M2—1)-4af. (43) 


The methods of § 2 can now be applied to prove this result rigorously, and 
also to show that the pressure on the cone, as calculated from the linearized 
values of f and f’ by equation (33), without neglecting the term f”, is 
asymptotically correct as a > 0. This work is here omitted. 

The angle between the conical shock-wave and the ‘Mach cone’ to which 
it tends as a > 0 is 


tan-!n—tan-18 ~ (n—Bf)/(1+f?) ~ 3(y+1)?Me(M2—1)-4o4. (44) 


A comparison between this result and figures given by Taylor and 
Maccoll (2) for « = tan 10° = 0-176 is adjoined. High accuracy is 
not to be expected (i) since a is too large for a factor 1+ O(a) to be 
negligible, (ii) since the numerical results are simply differences of numbers 
given to one decimal place of a degree and so may be considerably in 
error (thus the 0-6° value seems much too small according to differencing), 
(iii) since the form of (43) indicates that its accuracy falls to nothing as 
M tends to 1 or ©, a fact borne out by the rigorous treatment. It is 
therefore not believed that the discrepancies shown in the table are an 
index of any error in the work. 





M T.and M.\ Eg. (44) 


1°07 ize) 4° 
1°22 06 12 
1°39 08 I‘o 
1°81 o'9 I*2 
2°39 I°5 2°2 
3°33 2°6 51 








The strength of the shock-wave is usually measured by the relative | 
pressure oss (p.>—p,)/p,. This is found to be asymptotically | 
3(y+1)M*%(M?—1)-1a4. The deflexion produced by the shock-wave is | 
asymptotically mall )M4(M?—1)-!o4, an interestingly small amount 
compared with the total deflexion (~ «) of a stream-line near the cone. 
The entropy change at the shock is proportional to M8(M?—1)*s*. 
which for medium-sized supersonic M and small « is very small indeed. | 
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The most important conclusion from the results of this section can now 
be drawn. When air flows past any body of revolution whose tip is conical 
but which curves back convexly behind the tip, the shock-wave will, near 
the tip, have the same strength as if the whole body were conical with the 
same semi-angle, but farther out will not be so strong since less work need 
be done on the air farther out. Hence the maximum strength of the 
shock-wave ahead of any such body of revolution is known. In particular 
the maximum change of entropy is known, and, when this is very small 
indeed, the use of the adiabatic equations behind the shock-wave is fully 
justified. 


4. The uniform spherical wave 
When spherical instead of cylindrical symmetry is assumed, the only 
alteration in the equations and boundary conditions of § 2 is due to the 
change in V*¢. This doubles the term 2-1f’' in (7) and so doubles the term 
dz/dy in (20), and (25) then becomes 
y 
zy gi(y) = ygi(y)h(y) | ygi(y)h'(y) dy+ constant 


x 


y—93(y)h(y)— (log y)gi(y)h' (y) + 


| (log y)g}(y)[h"(y)+49(y)h'(y)| dy+K; (45) 


¥ 


(26) becomes z = (1+ Ky)|1+O(ylogy-)], so that the boundary con- 
dition z = a~®, y = a gives K ~ a. The estimation of the integral when 


y = € follows a similar course to that in § 2; and (29) then becomes 
M-* = h(e)—(e log e)h'(e) + O[(M—1)log( M—1)-*]+ 
+ O| (M— 1)log a-*]+- Of (M— 1)a“log(.M—1)-*]+ 
LK .4(y+1)(M—1)[1+0(M—1)], (46) 
or (using (22) and (24) to estimate h(e) and h’(e)) 
4(M—1)|log }(M—1)-*+-$][14+ 0(M—1)+ O(a4)] 
- 4(y+1)2K(M—1)1+0(M—1)+O(loga—)], (47) 
or, finally, log(M—1)~ —- - ~— —=» (48) 
y+1 (y+1)a8 
as a > 0. 
Thus the strength of the shock-wave is of an order of smallness 
exp|—(y+1)-1a-8]} rarely encountered in physical problems. This is 
borne out by Taylor (3) who gives a value ‘of order 10-1*’ when a = 0-2 for 


the strength (which is closely proportional to M—1); now log(10-19) = —44 
and —(y+1)-10-2-3 52, a close enough correspondence. Possibly 
5092.3 


Y 
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a more accurate approximation to M—1 could be found than a mere 
asymptotic value for its logarithm. It is so small, however (for the values 
of « which justify the use of approximate theories), as to be negligible in 
practice, lending little point to such an investigation. 

The asymptotic correctness of f on the sphere as deduced from linear 
theory, and hence of the pressure on the sphere as calculated from this 
by the full Bernoulli equation, now appears as follows: 

a a 
f(x) = | ydx = [xy],,— | x dy 
M € 


x 


= a+ O(e)—[1+ O(alog a-1)| [ (l+a-y)-? dy = —a?+ O(a8log a7), 

. (49) 

In the general case of waves produced by the expansion and contraction 

of a small sphere we may assume that the same result holds, i.e. that, if 

the linearized form of the potential 6 = r—!f(r—at) is used, then the more 
exact form p = p,—p(é¢/et+4 9?) of Bernoulli’s equation, i.e. 


P—Po = plarf'—i(r-4f'—1“*f P+ arf I}, (50) 


must be employed if pressures near the sphere are to be obtained correctly. 
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U+2% 


SUMMARY 


— 
~ . € . 7 . 
The function f(x) = du is tabulated to four decimals for the range of 
¢ &swzr 
rgument 2x 0(0-02)2(0-05)3(0-1)10. To facilitate interpolation for small x2 an 
ixiliary function g(x) f(x)+log x is tabulated to four decimals for x = 0(0-1)1. 


\n asymptotic expression is given for values of x greater than 10. Full details 
fthe method of computation are given and an interesting application of the Euler 


transformation to the summation of asymptotic series is included. 


|, Introduction 
Iya paper being prepared for publication, entitled ‘The Rectification and 
Observation of Signals in the Presence of Noise’, R. E. Burgesst has 
msidered, among other similar problems, the determination of the 
response of a detector to a random noise voltage having a narrow spec- 
ium. The relation between the output voltage V and the amplitude A 

{an applied signal is assumed to be of the form 

VocA for large A, 

x A* for small A, 


t,generally, V = A?/(A+C), where C is a constant. Various quantities 
lepresenting the response of such a detector to random noise are expres- 
‘ible in terms of the error function and other tabulated functions, and in 
ddition, the integral 


eo 


a 2 
e-u 


f(x) = | ——du. 
ute 
0 
This integral is of sufficiently simple form to be of general interest and 
twas thought worth while, when carrying out the more limited computa- 
ons required for the application mentioned above, to compute a more 


tailed table in which, for example, interpolation would be linear. This 
lable, together with a brief account of the methods of computation 


¢mployed, is given in this paper. 


' We are indebted to Mr. Burgess for the information contained in this paragraph. 
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2. Range and accuracy of the table 


The integral f(a) is tabulated to four decimals for the range of argument 

x = 0(0-02)2(0-05)3(0-1)10. A table of the auxiliary function 

g(x) = f(x)+logx (1 
is included to facilitate interpolation near x = 0, where the integral 
diverges. g(x) is tabulated for the range x = 0(0-01)1. 

The error in the last figure tabulated is less than 0-6 of a unit. For 
values of x greater than unity the table of f(x) is linear; that is, an inter- 
polated value when rounded off does not differ from the correct rounded- 
off value by more than one unit. For values of x less than unity it is 
intended that the table of g(x) should be used for interpolation and this 
table is linear except in the single interval 0 < # < 0-01. In this interval 
the error of an interpolated value can be two units, but this is quite 
negligible when account is taken of the infinity in f(7) when 2 = 0. 

It may be noted that the interval 0-00-0-01 is not the only one in which 
the second difference 5’g exceeds 4 units, which is generally taken as the 
upper limit for a linear table. That the table is nevertheless linear is due 
to the favourable nature of the rounding off in the pivotal values. 

For values of x greater than 10, f(x) is given to four-decimal accuracy 
by the expression 


; \n(1 1 1/1] ; 


3. Method of computation 

Different methods of computation were used for different ranges of the 
argument. For x less than 1 ascending series were summed, between | and 
4 a differential equation satisfied by f(x) was integrated numerically, and 
for x greater than 4 an asymptotic series was used. 

3.1. The asymptotic series 

If the denominator of the integrand is expanded in descending powers 
of x, f(x) can be written in the form 


fo) ioe) 








n—1 - 
r n n 
—_— 2 =— Uu 2 
f(x) = > (—)' ure“ du+ (—) | Se du 
ert an Uutx 
r=0 0 0 
n—1 
1 —) W/r+l1)\ , 
=;> (=! p(t) Re, 
2 grt 2 
r=0 
a 
1 £ . 1 n 
where |\R) <— | w*—te-" du = r{—}. 
x 2a" \2 
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| The series 


sument tony — XM(L 1 1-3 1.3.5 ta (3) 
J\ 2 \> 2x3 475 Sa? a DNa2 | at | ob! 78 | ty 


is accordingly asymptotic. (The first four terms of this series are quoted 





intecra] | it the formula (2) as giving four-decimal accuracy in f(x) for x > 10.) 
3.2. The differential equation 
it. For} From the identity 
n inter- ai 
- ue ] lf “ 
yunded- du scl : du 
' La» Ir 9 1 @»\8 
itv it is J u+a ae UZ J (u-+2) 
i 0 0 
ind this], a:¢¢ . . 
| the differential equation 
interva If 1 
( . 
is quite —— 2af N77 —— (4) 
dx 4 
| s obtained. 
n whiel The series solution of this equation in descending powers of x is, of 
as the . : , . . 
1 as the | course, the asymptotic expansion (3) already determined. Ascending 
r 1s due | series solutions can also be obtained, but for this purpose an examination 
of the behaviour of f(a) as x > 0 is required. 
ecuracy aa 
; 3.3. Limiting form near x 0. 


[It is apparent from the differential equation (4) that 
f(x)~ —logx as x>0 
and the limit of f(x)+logx, that is, of g(x), as 2 > 0 is of importance in 
formulating the ascending series. One method of determining this limit 


s of the | is as follows: 


nl and . a 
| du _— log t+ hax 
y, and 3 = ; 
. J (u-+1)(u+2) 1+<2? 
0 
y x 
. u oe ° 
a ( log «—4za . ] du 
[hus -du + — °- = | ew — —__ ]-—— 5 
powers J uta 1+-a2* 3 u*+-l/ut+a 
0 0 





and so, letting a > 0. 
, 2 1 \du 
lim{ f(a)+log ax} = | ad 
be ; 7 ( u®+1} u 
0 
] [ ‘ 1 \dt 
= = On ee ee FI oe 
2 | t+1/ t 
0 
aalie . 
3Y> (5) 


the integral for Euler’s constant being a well-known one due to Dirichlet. 
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3.4. Ascending series 





The singularity at « = 0 can be removed from the differential equation 
(4) by writing it in the form 
d 


2 , r2 ] r2 
—— ftlogx) = vre << —1). 


Integrating and using (5), an ascending series is obtained in the form 


= -2n+1 = -2n 7 
fte“loga = | vz > celine > ae by |e-z (6) 
_ n! (2n+1) —ni2n “ 


0 





An alternative form can be obtained by solving in series the differential 
equation satisfied by y = f+e-“logx. This equation is 


] s 
dy +2xry = vx—-(l—e-*) 
x 


dx 
’ (—)"a2"-1 
= fies ———= = ™ 
n! 
1 
@ 
The series solution y = >} a,x” has, from (5), a leading term a) = —}y 
0 
and it is readily obtained as 
@ = a) 
; ‘ l __\n (mn een : inl Ny2nt+1 a 
fte-@loga = = > Se is, 
9 ' i 6 - 5) ' 
o & n! ; 1.3.5...(2n+1) 


, : . @ 
where (2) is the digamma function Tn log T(x). 
da 


3.5. Details of computation 

The auxiliary function g(x), easily obtainable from f(x)-+e- log a, was 
computed from the ascending series (6) or (7) (occasionally both series 
were used for checking purposes). For the range 0(0-01)0-2 each individual 
value of g(x) was computed to six decimals. Between 0-2 and 1-0 six- 
decimal values were obtained at an interval of 0-05, and subsequently 
subtabulated to fifths on the National Accounting Machine. All values 
were then rounded off to four decimals. 

The main table of f(x) was computed in three stages. For x < 1, f(z) 
was obtained trivially from g(x). For the range x = 1(0-1)2(0-2)4 six- 
decimal values were computed by numerical integration of the differential 
equation (4), checks being applied at intervals by use.of the asymptotic 
series. These values were then subtabulated to the final intervals, chosen 
so that the table was linear above x = 1. For 4 <x < 10 the asymptotic 
series was evaluated to six decimals at an interval 0-5, and subtabulated 
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to fifths. Finally all values were rounded off to four decimals. For x > 10, 
the expression (2) is sufficient to give four-decimal accuracy. 

The method adopted for summing the asymptotic series for values of 
x less than 4 is of interest. The application of the Euler transformation+ 


y (—)"a S aw ed A"a 
ho — On+1 r 
oe n=0 ~ 

to the summation of convergent alternating series is well known; it is 
not so well known that the same transformation can often improve the 
accuracy obtainable from asymptotic series in which the terms alternate 
in sign. The accuracy obtained by summing an asymptotic series is of 
the order of the first term neglected. If this term is the smallest in the 
series, the succeeding terms form a divergent series which may be regarded 
summed’ by applying 
the Euler transformation, the transformed series being sometimes con- 


as representing the remainder. This series may be ‘ 
vergent and sometimes asymptotic. In the latter case the accuracy 
attainable may be improved by a further ‘Eulering’ and so on. 

As an illustration, consider the present asymptotic series for x = 1 for 
which the ascending series gives f(1) = 0-60513365. The third term is 
the smallest, with a value +0-4431135, so that, although in this case 
it has been shown that the remainder is actually less than the first term 
neglected, direct summation of the asymptotic series cannot give a better 
result than 


f(1) = 0°3862269+ R, where R < 0-44.... 


If the Euler transformation is applied to the divergent series of terms 
representing R, however, a new series, shown in Table 1, is obtained. 

It will be seen that the new series has a smallest term, the twelfth. 
Summing the first eleven terms, a new estimate for R is obtained as 


R = 0-2192169-+ R,, 


where R, is expected to be of the order of —0-0006. In fact this estimate 
for R gives f(1) = 0-6054438, which is correct to three decimal places, 
a substantial improvement. Still further improvement can be obtained 
by ‘Eulering’ the divergent series for R,, for which the differencing is 
shown in Table 1 and the traiustormed series in Table 2. 

To seven decimals the seventh term of this series vanishes so that a 
confident estimate for R, can be given as the sum of the quantities in 
Table 2. This value is —0-0003102, giving the third approximation 
(1) = 0-6051336, which is actually correct to seven decimal places! 


+ See, for example, Bromwich, An Introduction to the Theory of Infinite Series, 2nd ed., 
62. 
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TABLE 1 


2155 

- 142216 
+ 134730 
cae 39297 
a 29003 
- 15907 
f 12044 
— 8950 
bi 7554 
_ 6666 
rs 6307 Al A2 A’ At AS As 
= 6274 

—29gI 
4. (—) 6565 — 326 

617 —67 
- 7182 393 —48 

IOIO 115 —34 
+ (-) 8192 508 82 +7 

1518 197 —27 
= 9710 705 —109 

2223 — 306 
+ (—) 11933 IOII 

— 3234 

— 15167 


With alternative and much less laborious methods available, this method 
would not, of course, be used to obtain f(1). It is given here in full merely 
to illustrate the great power of the method, which is particularly useful 
in extending the range of applicability of an asymptotic series when this 
does not overlap the range of convergence of the corresponding develop- 
ment of the required function in an ascending series. 


The work described above has been carried out as part of the research 
programme of the National Physical Laboratory, and this paper is pub- 
lished by permission of the Director of the Laboratory. 
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THE FLAT DELTA WING AT INCIDENCE, AT 
SUPERSONIC SPEEDS, WHEN THE LEADING EDGES 
LIE OUTSIDE THE MACH CONE OF THE VERTEX 


By GWENDOLEN M. ROPER 
(Imperial College, London, S.W.7) 
[Received 5 December 1947] 


SUMMARY 

H. J. Stewart’s (1) analysis using the linearized equation of supersonic flows for 
finding the lift of a flat delta wing at supersonic speeds when the leading edges of 
the wing lie within the Mach cone of the vertex, is extended to the case when the 
leading edges of the wing lie outside the Mach cone of the vertex. 

From the condition specified on the wing, Fourier expansions are found for the 
downward vertical induced velocity component w, outside and on the Mach cone 
of the vertex. The constant velocities on the wing, outside the Mach cone of the 
vertex, are deduced from the shock-wave relations across the Mach planes of the 
leading edges. Fourier expansions are found for the axial and transverse disturbance 
velocities u, v, the constancy of these velocities, in regions outside the Mach cone, 
being deduced from the vorticity relations. 

The methods of the complex variable, used by Stewart, are used to obtain the 
induced velocity components inside the Mach cone of the vertex. The solution for 
the velocity component w is found from the conditions to be satisfied on the wing 
and on the Mach cone of the vertex. Hence, by using the vorticity relations and 
the Cauchy-Riemann relations, the solutions for the velocities u, v are found. 

The pressure distribution on the wing is deduced, and the total lift calculated, 
by integrating the pressure difference over the surface of the wing. The lift coefficient, 
based on area, is found to be independent of the angle of the wing. The corresponding 


drag coefficient is also found. 

1. Introduction 

Ustne the linearized equations of supersonic flow, H. J. Stewart has (1) 
found the lift of a flat delta wing at supersonic speeds, when the leading 
edges of the wing lie within the Mach cone of the vertex. 

The present paper extends Stewart’s analysis to the case when the 
leading edges of the wing lie outside the Mach cone of the vertex. The 
methods used to obtain the induced velocity components inside the Mach 
cone of the vertex are specifically the methods of the complex variable 
used by Stewart. 

Fourier expansions are found for the induced velocities outside, and on, 
the Mach cone of the vertex. The constant velocities on the wing outside 
the Mach cone are found from the shock-wave relations. 

The resuits obtained are not new, but it is sometimes useful to use 
different methods to obtain old results, since these methods may be 
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applicable to the solution of other problems. This is especially the case 
when, as here, the forceful methods of the theory of functions of a complex 
variable are used. 


2. Formulation of the problem 
The z-axis is taken through the vertex of the wing parallel to 
the undisturbed stream, and axes of 
y and z complete a right-handed 
>y orthogonal system, the y-axis being 
in the plane of the wing. The speed 
of the undisturbed stream is U, its 
Mach number M, and the Mach angle 








\ ° 
\ p. = cosec-1M. The semi-angle of the 
dé 4 ! SS _- wing is w, and the (small) angle of 


a incidence is «; in applying the boun- 
dary conditions, it will be sufficiently 
Fic. 1. The dotted lines represent the 
Mach cone of the vertex. The Mach angle : é 

is 4; the semi-angle of the wing is w,. Ponents of the disturbance in the 


accurate to put z = 0. Velocity com- 


x, y, and z directions respectively are 
denoted by w, v, w. The thickness of the wing is ignored. 
Using the linear perturbation theory, it is easy to show that each of 
these components satisfies the Prandtl-Glauert equation: 


22p 22D 72 
a ee tM (2.1) 


a.> 


6x2 By? -* az?’ 
where P denotes the dependent variable of the equation. 
The velocity components are related by the irrotational relations 


ow ov ou ow ov ou (2.2 


oy oz’ Oz dar’ ox ey 
The condition to be satisfied on the wing, that is, on z= 0, x >0, 
ly| <atanwo, is that tana = —w/(U-+u), or, to our order of approxi- 


mation, a ; (2.3) 


In what follows, we shall use 8 to denote ,/(/?—1). 


3. Transformation of equations 

We shall first refer briefly to the Hayes-Stewart change of variable, 
for the region inside the Mach cone of the vertex. 

By using coordinates formally analogous to spherical polar coordinates, 
H. J. Stewart (1) transforms the Prandtl-Glauert equation (2.1) into an 
equation in the form of Laplace’s equation in spherical polar coordinates. 
For the region inside the Mach cone of the vertex, P is taken to represent 
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the disturbance velocity w, and the equation is further transformed into 


(6S) eP 0. (3.1) 


Cs 06? 


the form of Laplace’s equation in two-dimensional polar coordinates (s, @) ; 


1' 2(y?-++--2?) ]-4 
H tan-!(y/z), § P ’ Mm ja—® _y 5 ) *". tan 
J \u+1 x 


here 


We shall obtain the equation to be satisfied by the induced velocity com- 

ponents inside the Mach cone of the vertex in the form of Laplace’s 

equation in two-dimensional Cartesian coordinates. The transformation 
of equation (2.1) is made by putting 

x 3rcothe, y—=rcosechocos#, 2z==—rcosechosin#é. (3.3) 

It can be shown that if (2,, 22,23), (Y1, Ye, Y3) are two sets of variables, 


suc t] at 9 9 9 9 9 9 9 9 9 
ach tha dx? —dx3—daz = h? dy? —h§ dy3—hj dy}. 


ind if P is an arbitrary scalar function, then 


2P @P @P l (hah, @P" C h,h, eP C h,h, eP 
re ey! oxe h,hghs OY; \ h, ms OYs h, OYs  €Ys hs OYs ‘ 


(3.4) 


It is easy to show that 


) 


2 


(F) dy? —dz* dr? —r? cosech?a(do?+d6é?), 


ind therefore 

Wor oP @eP l Sis , oP 0 [eP @ (eP 

— : - - - x2 cosech*o — | —— | — | — —|—]} |, 
ox* cy~ Cz* r2 cosech?2c | er cr co\ co o6 ral] 


(3.5) 
x? 1 
r ,—(y?+ =|! 
E 


B(y?+-2° 


where. by (3.3). 


(3.6) 


i 
sech o = Btanw 


A 
6 = tan-1(z/y) 
v, 8 being the normal spherical polar angular coordinates. 

In a cone-field, the conditions are constant along any radius vector, so 
that the disturbance velocities depend only on non-dimensional coordi- 
nates, and the velocity components wu, v, w are independent of r. 

Therefore, when P represents one of the velocity components in a cone- 
feld, equation (2.1) becomes 

AP FP 
éo2 | 0b? a 


71s real within the Mach cone of the vertex, and is zero on the Mach cone. 


0: (3.7) 
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On the axis, y = 0, z = 0, r= 2/B, co = ©. 
For the region outside the Mach cone of the vertex, replace o by —iy, 


and use 
x = Breot x, y = rcosec x cos 8, z=recosecysin@; (3.8) 
we find 
pAP_AP_AP 1 [a(scoeoty®P)_ AP , API 
da® = dy* 2% ~=r* cosec*y | ar or] ex? oF 
—— (3.9) 
where r = [y2-+22—22/B}t | 
y? +22) 
sec x = Bly" Es Btanw }. (3.10) 
6 = tan-!(z/y) | 


Therefore, when P represents any velocity component in a cone-field, 
equation (2.1) becomes 

52P e2P 

é C 

eon sean sm WE (3.11) 


9 


6x? oP? 

x is real outside the Mach cone of the vertex, and is zero on the Mach cone. 
On the leading edge, 

sec x = sec xy = Ptanwy. (3.12) 


Solutions of the equations for the regions outside and inside the Mach 
cone of the vertex will now be obtained. 


4. Solution for the region outside the Mach cone of the vertex 
The equation to be satisfied by the velocity components wu, v, w for the 
region outside the Mach cone of the vertex is 


i (3.11) 
ox? 0? 
of which the general solution is 
P = f(6+ x)+F(0—x). (4.1) 
The condition to be satisfied by w is 
w= Wy = —Ua (4.2) 
on the wing, i.e. for @ = 0 or 7, 0 < |x| < sec—(Btan w,). 


There is symmetry about the x-axis, and the value of w, as a function 
of 0, should be the same as the value of w as a function of ~—6. Also the 
velocity w has the same sign on the two sides of the wing. 
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Therefore, for the region outside the Mach cone of the vertex, we assume 


the Fourier expansion 


w = Ay+} ¥ [A,,{cos 2n(8+ x)-+-cos 2n(0—x)}+ 


+ Aj, {sin 2n(@+ x)—sin 2n(0—x)}] 


x 
Ay+ > [(Ag, cos 2nx+A3, sin 2nx)cos 2n6]. (4.3) 
n=1 
On the wing, @ 0 or 7, and 
a 2 
w= Agt > Ap, cos2ny+ > Az, sin 2nx. (4.4) 
n=1 n=1 
The condition (4.2), to be satisfied by w, gives 
U xX Xo w Wy = — Ua 
| 6=0or7. (4.5) 
Xo X har, Ww 0. | 
Therefore me ow 
= dy <9 Xo 
A, = | wdyx E 
0 
4 2w, sin 2n 
ro w cos 2nx dx - ae ee Xo. 
. w 7 n 
0 
— | 2w,[1—cos 2n 
A,, wsin 2nx dx —_— — 0 Xo ; 
aa” 7 n 
0 
On the wing, 
Ua — /sin 2nx , sin 2n(x»—x) 
Hecate ee NK & . 
u ye xo ( n — + - a ’ (4.6) 


n=1 


nd for the region outside the Mach cone of the vertex, 


T 
n n 


4 


2U 4 , xo [sin 2nyxcos2n6 | sin 2n(yy—x)cos 2n6 7 

—"|xet > ( + ey. (47) 
n 1 

Un the Mach cone of the vertex, x = 0, and 


2U a S (= 2x9 COS my 
n : 


w | xo+ 
n=1 


summation the 


(4.8) 


7i 


Abel’s method of Fourier series in 


Summing by 


equation (4.7), we obtain, for the region outside the Mach cone of the 


Ua ) Li] ] —e2(8+(xo-x)}4 
7 2(xo—x) +1 log | —eXO—(xo-x)}4] |’ 


except at |@ Xo—xX; Or |0 


vertex, 
(4.9) 


w Re- 


™—(xo—x), Where there are singularities. 
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Taking the principal values of the logarithms, we obtain the values 


of w: 
; . __ —Uar, ” 
0 < |A| < xXo—x;, asian —[2(xo—x)—2(xo—x) +7], 
™—(xXo—x) < |0| <2, - —Ua. 
—Uea ‘ 
Xo—X < |! < w—(xXo—x); lias —[2(xo—x)—2(xo—x)] = ¥, 
(4.10) 
The values of w on the Mach cone of the vertex are given by 
o= 10] =<», } 
io ee 
a — Ke * Ol<a J (4.11) 


Xo - 6 = Mas w= U0. 


The surfaces xy = constant are cones, coaxial with the Mach cone of the 
vertex, and the surfaces 6 = constant are axial planes. The surface y = 0 
is the Mach cone of the vertex, and the surface x = xp is the cone through 
the leading edges of the wing. The surfaces 0 = + x9, 6 = +(7—xp) are 
the axial planes normal to the Mach planes of the leading edges (see Fig. 2). 
Therefore, we have deduced that, in the regions bounded by the Mach 
cone of the vertex (|! < x9; 7—x9 < |0| <7) and the Mach planes of 
the leading edges, w is constant, and equal to —U«, that is, equal to the 
value of w on the wing, and that everywhere else outside the Mach cone 
of the vertex, w = 0. 

Cw Ov ow 


Mm = , cou ; 
rhe vorticity relations — = — ,— = — give 
Cz Ox Oz oy 


in y sing l/ow cw 
8 J —S = — | 
Oy os Blog op 
sin — — cos po” — cos” 

us ' @ é 


\ (4.12) 
° 7 Ov 
—sin UY) | 
where 4 = 0+y, fb = 0—x. 
Writing the solutions of equation (4.12) for the velocities uw, v in the 
forms 6 = wld)-+- wld), 
v = o(f)+v(), 


since the relations (4.12) are true for all values of % and 4%, 


sing 2?) - 1 ow 


(4.13) 





| Oh Bop (4.14) 
sin oe = cosy 
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nd there are similar relations for u(y) and v(%). These relations are 
snalogous to the vorticity relations (5.7, 5.8), used later for the solution 








Fic. 2. A4OB isa Mach plane of the leading edge OA of the wing. From the geometry 
fthe figure, the angle subtended at the axis, OC, of the Mach cone of the vertex 
y the Mach plane AOB is given by 1/8 = tanp = tanw,cosACB. Therefore 
ingle ACB = yp, since B tan wy = sec xp. Angle ACD = yo—x. 


We have shown that w = 0 or —Ua everywhere outside the Mach cone 
f the vertex; therefore 
ow ow 
ous . orf i 


everywhere, except at the singularities at 0 = +x 5, 0 = +(7—xo), where 


0 


the derivatives are undefined. Therefore, from relations (4.14), we should 
expect wu and v to be constant in the regions outside the Mach cone of the 
vertex, except perhaps at the singularities. 

CW OW ° . 

ab Bd are zero everywhere outside the Mach cone of the vertex, and 
the series obtained by differentiating the Fourier expansion for w, term 
ty term, do not converge. But it can be verified that the Fourier expan- 
ions for w« and v, which we shall deduce by using physical arguments, 
utisfy the relations (4.14) when they and the series for w are differentiated 
term by term. By use of a convergence factor, a rigorous demonstration 
8 accessible. 

The values of the constant velocities on the wing, outside the Mach 
one of the vertex, can be deduced from the shock-wave conditions across 
the Mach planes of the leading edges. Assuming the continuity of the 
velocities up to the wing, in the regions between the Mach cone of the 
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vertex and the Mach planes of the leading edges, the values of these 
velocities will also be known. Outside these regions, wu, v, w are zero. 
The shock-wave conditions, across the Mach planes of the leading 
edges, are that the tangential velocities are continuous. 
Therefore for 0 < 6 < x9 (u = Up, Vv = Vp on the wing for 6 = +0) 


UB = (U+tug)B-+vy 00s xo-+ty sin Xo | 


(4.15) 
U cos wy = (U+u9)c0s wy+ vy Sin wy J 
Therefore, remembering that sec yy = B tan wy, 
—W, Ua 
uy = = 
Bsinx,  fPsiny, }. (4.16) 
y= U x cot Xo 


Similarly it can be shown that, on the wing, 





for 6 = —0, % = —e ; v = Uacot x, 
Bsin xo 

for 6 = —7, “= - es v = —Uacot xp. 
Asin xo 


Therefore the conditions to be satisfied by u, v for the region outside, and 
on, the Mach cone of the vertex, are: 





0 <6 < xo—x: “= Me, Y == Ue, 
7™—(xX9—X) <9 <7, © == We, v= —UVp, 
—(xo—x) < 8 <0, u = —Up, v= —Up, } (4.17) 
—7 <0 < —[7—(x9—x)], u = —Up, VU = %, 
Xo—X < |8! <[a7—(x9—x)], «= @, v= 0, 
where uy) = Ua/(Bsin x9), vy = —Uacot xp. 


Assume the Fourier expansions: 


a 


u= > [B,,.,sin(2n+1)6], 
n=0 

v= y [C,,, sin 26], 
n=1 


where B , C,,, are functions of y. Hence, we have 
2n+1 2n X 





- 4Ua [1—cos(2n+1)(yo—x 
Bo4s i | usin(2n+1)6 dé = - 3 Ua a Mixon ; 
J 7B SIN Xo 2n-+ 
, _ —2U acot xq j= =H 0) 
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| where U, V, W are analytic functions of the complex variable ¢ = se‘. 
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U . —cos{(2n+ a . 
= fe > [font Dew \sin(2n+1)8], (4.18) 
mB sin Xo — 2n+1 J 
F -2U a cot x ys |= ike sin 2no| (4.19) 
- Zoe a8 tt , 
n=1 


It is seen that w, v are of the form given in (4.13), since the Abel sums of 


SS sin(2n+1)0 S sin 2n0 


/ 


haat 2n+1 n 
n=0 n=1 
are respectively equal to +}7, +($7—8@) = +(5—*"), according as 
his2 0. 
On the Mach cone of the vertex, y = 0, and 
| 4Ua << [{1—cos(2n+1)xo) .. 
u <i. i, * js — oom )Xo\ sin(2n+-1)8 . (4.20) 
7B sin Xo My 2n+1 | 
. 2U «cot xo +. {1—cos 2X0) oi ond |. (4.21) 
7 hat \ n | 
n=1 
There are singularities at 0 = +(x ,—x) and @ = +{7—(x9— x)}. There 


are also discontinuities in the values of uw, v at @= 0 and 6 = 7, that 
is, at the wing, where the sign changes. 


The values of u, v, w on the Mach cone of the vertex are given by: 


0<0 <x. u Ua/Bsinyy, v= —Uacotyxy, w = —Ua, 
TX) <A <7, u = Ua/Bsinyy, v= Uacotx, w= —Ua, 
—x¥, <9 <9, u = —Ua/Bsin x5, v = Uxcoty,, w= —Ua, 
—7 <6- (7—Xp), U Ua/Bsin xo, v = —Uacot x9, w = —Ua, 
Xo * 6 T™— Xo; u 0, vy = Q, w= 0. 


(4.22) 
It will be seen later that these values of u, v on the Mach cone of the 
vertex agree with those obtained from the solution for the region inside 


the Mach cone. 


5. Solution for the region inside the Mach cone of the vertex 
when the leading edges of the wing lie outside the Mach cone 
For the region inside the Mach cone of the vertex, H. J. Stewart finds 

the solutions to equation (3.1), for the case in which the leading edges 

of the delta wing lie inside the Mach cone. He uses the complex functions 


U = uid, V= 0+ 78, W = wie, (5.1) 
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By using the Cauchy-Riemann relations for P, whose real part is w, v, or w, 
and the vorticity relations (2.2), he obtains the relations between U, V, W: 
dV ii — " dW 


dt — 142 dl - 
dU aw sa 
dg B(1+- 2) d& 


The solution for the region inside the Mach cone of the vertex, when 
the leading edges of the wing lie outside the Mach cone, can be obtained 
equation (3.1) for this case, or by solving equation (3.7) in terms of the 
complex variable # = 6+-i0. 

The solution is obtained here from equation (3.7), this equation being 
analogous to equation (3.11), used for the solution for the region outside 
the Mach cone of the vertex. 

Since the velocity components, wu, v, w satisfy the equation 


or oP 
-+—— = 0, (3.7) 
Go? 06? 
the complex functions U = u+7a, V = v+i0, W = w+ia, are taken as 


analytic functions of the complex variable # = 6+ic. The Cauchy- 
Riemann relations for P are 


oP @aP oP _ oP 53 
—_ = —, ——. 5.3) 
06 oo Go ali] 
The vorticity relations (2.2) give 
ov - ,ow Cw i. wee 
cos 0— +-sin@— = cotho|cos#@— — sin #—}, 
06 06 06 06 
5 Oe ou ] Cw a 
coth o sin @— + cos 6— = —-—cosecho—, (5.4) 
00 o~—Sté«@B’ 06” 
ou i Oe ] oo 
coth o cos 6— — sin 6— cosech o- 
06 7 an  B 
The linearized equation of continuity 
ov cw ou 
oe pee 
Cy =z 
gives 
d ‘ee : Ow ou 
coth o cos 6 — sin @— + cothosin 6— oT “ = —Becosecho—. 
06 06 06 + 6 = . 06 
(5.5 


From equations (5.4) and (5.5), 


by extending Stewart’s use of the complex variable { to the solution of 


aD . ,OW noe ow ep 
cos @— +- sin @— = cotho|sin @— — cos —}. (5.0) 
06 06 06 06 
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The first of equat ions (5.4) added to 7 times equation (5.6) gives 


ev .  .e€W ; . »ov ow 
cos 6 | sin @ i coth of sin 6 — — cos 0 — }. 
ob ral] Paley fale 
Hence, since V and W are functions of #, 
dV i coth o cos 6+ sin 6 dW dW — 
_ = cot? —. (5.7) 
dd cos 6—icothosin# dé dd 
Similarly it can be shown that 
dU ] aW - 
; — = cosec § —, (5.8) 
dd B dd 
Inside the Mach cone of the vertex, 
0 o¢<c 
On the Mach cone, o = 0, ~~ <O<-«. 
On the wing, 6—O0orn. 
00 
AA AF 











- -(™-X5) —Xo 0 Tt TT 
(7X0) -X Xo Xo 


B C D —E 6 


Fic. 3. #-plane. Points A, B, C, D, E, F in Fig. 3 correspond 
to points A, B, C, D, E, F in Fig. 4. 


On the d-plane, the Mach cone is represented by the 0-axis from 0 = —7z 
to 6 +o, and the part of the wing inside the Mach cone by the two 
semi-infinite lines 6 = 0, 0 = 2 (o > 0). The transformation 


C = §+in = —cos 3 (5.9) 


transforms the semi-infinite region bounded by the lines 6 = 0, 0 = z, 


7=0(0 <6 < 27,0 > 0) in the #-plane, into the region 7 > 0 in the 
‘plane. The contour is mapped on the whole é-axis. 

The points —1, —cos xy», Cos xp, +1 on the €-axis correspond to the 
points 0, yx», 7 Xo, 7 on the 6-axis. 


Che conditions to be satisfied by w inside the Mach cone are: 


w w Ua on the wing, 
0 t =) 
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and 
I< |4| < : 
es Xo w= —Ua = wy 
™—Xy <|Ol| <a on the Mach cone. 
— | a — " 1= Q 
Xo <|8| <a7—xy, u (5.10) 
An 
— —1 —cos Xo | 0 COSX%q =f ~ s 
A BOC DE .* 





Fic. 4. ¢-plane. Points A, B, C, D, E, F in Fig. 4 correspond to points A, B, C, D, 
E, F in Fig. 3. 


Therefore, for the region considered here, the conditions to be satisfied 


are: 
—o <€ < —COS xy , 
w= —Ua : 
COS Xp < € < +00 on the é-axis. (5.11) 
—COS x9 < € < COS x9, w= 0 


There is discontinuity in the value of w when é = -+cos yy. Conditions 
(5.11) are satisfied by 
W = +=" Wot Iog{S—008 Xo| (5.12) 
©\C-+-e08 xo)’ 
the principal value of the eel being taken. 














Therefore W = Wy+ 0 Jo g/— er —-e (5.13) 
|—cos 3+ cos xy 
on the upper surface of the wing. 
IW = 2wyt = cos x,sind . 
Hence feae ie oe ee (5.14) 
dd 7 COS*d—cos?xo, 
Therefore, using relation (5.8), 
dU a __ 2wot COS Xo . (5.15) 
dd Br cos?d—cos*y, 
Wet {tan Xo—tan i) 4 (5.16) 


~ Brsiny, \tan xo-tand| - 

the principal value of the logarithm being taken , and C being a constant. 

Above the wing, 0 <6 <7, and the conditions to be satisfied by 4 
above the wing are: 


0<0< 
” u = constant = u,, : - 
T—Xy9 <O< 7 fora = 0. (5.17) 


Xo <9<7—-xy, w=O 
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For 0 <@ < Xo, or 7—Xy <8 <7, 
U —sinh 2oe \ 


xX 
“= scenes tan-1/ $$} 
Br sin xo | (cos 20+ cosh 2c)tany,—sin 26) 


oe sinh 2o 40 
\(cos 20+ cosh 2c)tan xy+sin 24 rr 
For xp < 9 < 7— xp, 








Ue al —sinh 2o | 
“u = ~——— | tan :; —__——, —7— 
Br sin Xo \(cos 20+-cosh 2c)tan x,—sin 26} 
eank, O 
-tan-1/_— _ sinh a . \ C. 
\(cos 26+ cosh 2c)tan y,+sin 26} 
U Uc 
Therefore C ons ee t, = C= a (5.18) 
Bsin xo SIN Xo 


On the upper surface of the wing, 6 = 0 or z. Therefore 


sae Ua —_ if sinh2o _ vw sinh2o | a* 
Br sin xo \tan y»(1-+cosh 20)} \tan x9(1-++cosh 2c)/ 
77 
Br sin Xo 
where sinA:cosA:1 = cot xy: coth o: +,/(cot?x + coth*e). 
9717 ¢ 
Therefore u * _ sit -1f - _tan Xo \ 


= —sin~*{—— AS. } . 
Br sin xo \,/(tan*x»+tanh?o)} 
or, putting ¢ = sechocos xy (¢ < 1), 
OTT 2 
— =. ee (5.19) 
Br sin xo (1—#?*) 


on the upper surface of the wing. 


Uu = 


On the under surface of the wing, 

OTT hd 
—2U0 ent, Xo 
Bar sin xo /(1—#*) 
Using the relation (5.7), 

dV 2wyi cos x,cos? 
dd TT sin?y,—sin??’ 
U cei sin #—sin 
Vv oct a eae in Xo\ (5.21) 
- \sin #+-sin Xo! 
The conditions to be satisfied by v above the wing, for 0 < 6 < 4z, are 
0<0< xX, =x%\v= constant = v,, 


for o = 0. (5.22) 
Xo * 0 < 4n, v Q 


Therefore D 0. v. = —Uacot Xo (5.23) 


m 
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Equation (5.21) then gives, for 7—y,» <9@<7,0= 0, 
v = +Uacot xp. 


On the upper surface of the wing, for 0 = 0, 


Ua sinho sinh o 
v= —— cot xy tan-1}- moe | e—tan-/ = | 
7 |—sin xo) | sin xo) 
2U« . 
= — = cot xo($a—2’), 
where sin X’: cos A’: 1 = sinho:sin x9: +,/(sinh?o+-sin? x). 


Putting t/ = cosechcsin xp, 


2U , 
gon _ 20 x cot | a. _ (5.24) 
7 \ (1+#*) 
On the upper surface of the wing, for 0 = 7, 
a Pe ua’, - 
T \ (1+#°) 


On the under surface of the wing, v has the same numerical values, 
but with a + sign for 6 = 0, and a — sign for 0 = =. 

The values obtained for u, v on the Mach cone (5.18, 5.23) agree with 
those obtained in (4.22). 


> 


im ure 


n B 0 Cc 





Fic. 5. ¢-plane. Angle COD = angle AOB = yo. Points A, B, C, D, E 
in Fig. 5, correspond to points A, B, C, D, E in Fig. 6. 


—see xy -I 0 +1 SOX coy! 


AB Cc De 


Fic. 6. ¢’-plane. Points A, B, C, D, E in Fig. 6, correspond to points 
A, B, C, D, E in Fig. 5. 





It may also be of interest to compare the solution of equation (3.1) with 
Stewart’s solution for the case when the leading edges are inside the Mach 
cone of the vertex. 





Us 


putt 
01 
cirel 
the 
n-ar 


T 


tran 
regi 
fa. 


T 


Fro 
Va 


He 


Th 
of 
vel 





ralues, 


» with 


) with 
Mach 











DELTA WING AT INCIDENCE, 





FLAT AT SUPERSONIC SPEEDS 341 
Using Stewart’s notation, and taking 6 = +47 on the wing, that is, 
putting tan@ = y/z, U, V, W are taken as analytic functions of £ = se‘. 
On the ¢-plane, the Mach cone of the vertex is represented by the unit 
circle with centre at the origin, and the part of the wing which lies inside 
the Mach cone, by the diameter of the unit circle which lies along the 
n-aXxis. 
The transformation 


e’ Lin! (5.25) 


| 
us 
nw 


transforms the interior of the unit circle, in the region € > 0, into the 
region 7’ > 0 on the ¢’-plane. The closed contour is mapped on the 
é'-axis. 


The value of W is given by 


iUa, {2t¢ cos yy—C?+1) Eas 
W — log; — = ‘ (5.26) 

T 2if cos yp + C?+1) 
From (5.26), and Stewart’s vorticity relations (5.2), the values of U and 


V are obtained, 


U (24 cos 2 
U a 2xo- 2tan-1/2 1008 Xo\ | (5.27) 
Br sin xo | sin2x, | 
V Wacol X0Jog/$—26sin Xor}), (5.28) 
T \€?+2¢sin yy +1) 


Hence the solutions (5.19, 5.24) for w and v follow, where in this solution, 


ds 


¢ — Sear 
ZS COS y P 28 S81n 
f Xo O a ae. 


i+? * 1—s* 
The value of dW/dZ when y, = 0 is in agreement with Stewart’s value 
of dW /d¢ when the leading edges of the wing lie in the Mach cone of the 
vertex. There is also agreement in the values of w for this case. 


6. The lift and drag of the delta wing 

Taking the disturbance velocities u, v, w to be small compared with the 
undisturbed stream velocity U, and retaining only their first powers, 
Bernoulli’s equation 


“dp P 
Ea iq? = constant, 
oe 
od ] 
pee ( ) y 
vives PsUu 0, 
J P 
Po 


Where p is the density of the undisturbed stream. 
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Therefore the pressure on a surface element of the upper surface of the 
wing is given by 


—2pU?a sin 
pe | OR selling” | cs 
p pUu — Eaves (6.1) 
inside the Mach cone of the vertex, and 
—pUx 
Ap = —f-— 6.2 
P Bsin xo Xo (6.2) 


outside the Mach cone of the vertex. 

The pressure on the under surface of the wing has the opposite sign. 
Therefore the pressure difference is 2|Ap]}. 

To find the total lift, the pressure difference is integrated over the 
surface of the wing. 








Fic. 7. OC = the maximum chord of the wing, c 
OR=R ai (a? + y?). 


On the wing, R? = 2?+ -y*, tanw = y/x, t = secha cos x, = B tanw os xp, 


and c is the maximum chord of the wing. 


Therefore the total lift 


csecw COS Xo 
4p U2x 2 
T 


Se sin-? ___*4. sin Xo dt+ és cos*wR dR 
B* sin x9 COS Xo 
R 


J(1—#?) 


=0 0 C08 x, 


n 

9 nTT2~¢2 9 r. 

= a = — | — [ sin Xo dt+-1—n|, 
B?n,/(1—n?) | 7 (1—#?) 











writing Cos x) = n. 
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n 
2 
Now [ PE oat a Xo — 
V(1—#) 2 (1—#),/(n?—#?) 
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= (n—1)4r+.,/(1—n?)}z. 
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Therefore total lift 


Io 2ac2 9 U2axc2 
=: a. . - = [n—1+./(1—n?)+1—n] = —. (6.3) 
Bn,/(1—n?) B*n 
The lift coefficient C;, based on area, is equal to 
the total lift dou (6.4) 


loU?.area of the wing f%cosx,tanw, 8 
Therefore C;, is independent of the angle of the wing. 
The pressure difference is finite at the leading edges. Therefore if the 
thickness of the wing is neglected, and we assume that the resultant force 
on the wing is normal to the wing, the total drag is 


2pU2a2c? 


arena 5 (6.5) 
B* cos xo 
The corresponding drag coefficient is 
‘ 4a? 
: ¢ 
ae B ° (6.6) 


7. Conclusion 

For the flat delta wing, at incidence, at supersonic speeds, when the 
leading edges lie outside the Mach cone of the vertex, the formulae for 
the disturbance velocities in the regions outside and inside the Mach cone 
of the vertex have been obtained. Hence the pressure distribution on the 
wing has been deduced and the total lift calculated. 

The formulae for the disturbance velocity u, and the pressure difference 
Ap, agree with those found, using different methods, by A. E. Puckett (2) 
and G. N. Ward (3). 

The lift coefficient C,, based on area, is found to be equal to 4a/8, and 
therefore independent of the angle of the wing, a result previously obtained 
by Puckett. 

The yawed delta wing, with the leading edges outside or inside the 
Mach cone of the vertex, will be considered in another paper. 
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included the reading of the paper, and the making of valuable criticisms 
and suggestions. 
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ON THE HODOGRAPH TRANSFORMATION 
FOR HIGH-SPEED FLOW 
I. A FLOW WITHOUT CIRCULATION 


By 8. GOLDSTEIN, M. J. LIGHTHILL, and J. W. CRAGGS 
(Department of Mathematics, The University, Manchester) 


[Received 13 December 1947] 


INTRODUCTION AND SUMMARY 

In ref. (4) a solution is given to the problem of finding, by use of the hodograph 
transformation, steady plane isentropic flows round contours, reducing to given 
incompressible flows when the Mach number at infinity tends to zero. The solution 
contains an arbitrary function and therefore gives an infinity of such flows for each 
incompressible flow, when circulation is absent: with circulation, however, only one 
determination of this arbitrary function is permissible. 

In this paper an alternative method of solution is indicated. This method is not 
so easy to put into a general form, so it is set out for the case when the incompressible 
flow is the symmetrical one about a circle. Various possible generalizations are 
indicated. The approach is presented as one which helps to clarify understanding 
of the use of the hodograph transformation for the study of the flow around bodies 
and as one which may suggest future investigations on the application of the trans- 
formation. It is shown also how the results obtained for the case studied follow 
from the method of ref. (4). (We understand that this case is among those studied 
by Cherry in a paper to appear in the Proceedings of the Royal Society;+ the method 
there used, however, is similar to that of ref. (4): our interest here is in describing 
the alternative method.) 

1. Statement of the problem 

In a steady, two-dimensional, isentropic, irrotational motion of a perfect 
gas with constant specific heats, with viscosity, heat conduction, and 
radiation, and the effect of body forces (such as gravity) on the gas 
neglected, it follows from the equation for the absence of vorticity and 
from the equation of continuity that there is a velocity potential ¢ and 
a stream-function % such that 


“= db, : Poy P; v dy, nn — pots p> (1) 
«x and y are rectangular Cartesian coordinates, wu and v the corresponding 
components of the velocity of the gas, and p the density of the gas; 
suffixes x and y denote partial derivatives and the suffix zero values at 
a stagnation point. 

If q is the magnitude of the resultant fluid velocity, c the local velocity 

of sound, y the ratio of the specific heats, and 
i— —., (2) 

y—l1 


+ This paper has now been published. See Ref. (6) at end. 
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HODOGRAPH TRANSFORMATION FOR HIGH-SPEED FLOW 
the equations of motion have an integral 
2Be?+-9g° = qq = 2Beq = (2B+1)q5, (3) 


where q,; is the (theoretical) greatest possible value of g, when c = 0, and 


q, is the value of g when q = c. 

With x and y as independent variables, it is possible to find a non-linear 
partial differential equation for ¢, or an even more complicated equation 
for ’. With uw and v as independent variables, however, there are linear 
equations for ¢ and ys; the equation for ¢% is the simpler. Instead of u 
and v, we use as independent variables 7 and @, where @ is the angle which 
the velocity vector makes with the axis of x, so that 


uUu—tv = ge-®, (4) 
and Tt = q?/q3.. (5) 


The equation for y is then 


472(] nes + 47/ 1+ (B— 1] + [18+ Ie] = 0, (6) 
which has solutions &, (r)e*'"’, where 
ys, (7) = 7?"F, (7), (7) 
F (7) is the hypergeometric function, defined by 
alt) = 1+ iat"? See - 
for 0 <7 <1, and 
a,+b, = n—B, a,b, = —4Bn(n+1), (9) 


2a, = n—B+[(2B+ 1)n?+B?}}, 2b, = n—B—[(2B+1)n?+B?}}. 
(10) 
If now in a low-speed potential flow past a cylindrical obstacle with 
undisturbed velocity U at infinity, we write 
z= x+iy, w= d+, C = qe-” — dw/dz, (11) 
then w is a function of ¢, and if we suppose that w may be expanded in 
a series w =>, (tjUy" (12) 
we may try to find a high-speed flow, ‘corresponding’ to the given low- 
speed flow, by taking the stream-function of the high-speed flow to be the 
imaginary part of a function W, where 
W = > d,f,(r1),(r)e-*", (13) 
7, being the value of r corresponding to the undisturbed velocity at infinity 
in the high-speed flow, and /,,(7,) being chosen so that (13) reduces to (12) 
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when gg > © and 7 and 7, > 0, ie. (since ¢,,(7) ~ 7?” as 7 > 0) so that 
f(t) ~ 77 as 7, > 0. 

An almost identical procedure was successfully used by Chaplygin (1), 
with f,,(7,) = [¥,,(7,)|“, to solve the problem of a subsonic jet of gas 
escaping through a slit in a plane wall. In the application to flow past 
an obstacle, however, certain difficulties are encountered. 

First, it is of most interest to find solutions in which there is a region 
of supersonic flow. If there is a solution of the form considered, with 
a single-valued function of 7 and @, branch lines (Lighthill (2) and Craggs 
(3)) cannot occur, but with a supersonic region present, limit lines (Craggs 
(3) and references there given) will appear above a certain speed, and the 
solution is then meaningless. 

Secondly, it is necessary to ensure that the streamline corresponding to 
the contour of the obstacle in the high-speed flow is closed. This can be 
done, but the shape of the contour varies as the Mach number of the 
undisturbed flow varies, and the problem of calculating the high-speed 
flow past a cylinder of given section has not been solved. 

Thirdly, even in the simplest case, w as a function of { for the low-speed 
flow has a branch point at the point in the ¢-plane corresponding to the 
velocity at infinity, and to represent w by series of the type (12) in the 
whole of the fluid requires four series. These series are analytic continua- 
tions of one another, but it does not follow that series formed from them 
in the same way as (13) is formed from (12) are analytic continuations 
of one another, and in fact they are not. If a series such as (13) is valid 
in any one region we are therefore faced with the problem of finding its 
continuations into the other regions of the gas. The simplest case seems 
to be that of the flow corresponding to low-speed flow past a circular 
cylinder, and it is the problem of finding series valid in the whole of the 
gas and analytic continuations of one another for this simplest case which 
forms the subject of this paper. 

When circulation is present the matter is more complicated still; w then 
has a pole at the point corresponding to the velocity at infinity, and a 
branch point elsewhere. 

Lighthill has lately obtained a very general solution, with considerable 
extensions, of the problem here considered. After considering the pro- 
perties of the functions %,,(7) (Lighthill, 5), he (Lighthill, 4) has shown 
how solutions for high-speed flows, valid throughout the whole of the 
region of subsonic flow, may be constructed in a great variety when 
circulation is absent, and also in a more restricted manner when circula- 
tion is present, independently of the existence of series such as (12) for 
the low-speed flow; these solutions may be transformed into series of the 
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type (13), and its correct continuations, when series such as (12) exist; 
and in the latter case the resulting series are valid also in the region of 
supersonic flow (in the absence of limit lines). We understand that the 
problem has also been considered by Cherry (6), in a paper which will be 
published in the Proceedings of the Royal Society. 

"Nevertheless it seems that publication of the solution of a definite 
problem (in fact the simplest possible) by a different and comparatively 
simple method is justified. The method used makes use of a contour 
integral of Barnes’s type to find the continuation of a given series, and 
isexactly analogous to the use of Barnes’s integral to find the continuation 
of the hypergeometric series (Whittaker and Watson, 7). The method 
makes it very clear that the reason why series formed in each region from 
the series for the low-speed flow in the same way as (13) is formed from 
12) are not analytic continuations of one another is that %,(7) has poles 
as a function of v; that the most convenient form for f,,(7,) is one for which 
f,(7,) as a function of v has no poles, and that, since y,(7) has zeros as a 
function of v, Chaplygin’s form for f,,(7,), although appropriate for the 
problem of a subsonic jet, is not suitable for the discussion of a high-speed 
flow corresponding to a given low-speed flow past a cylinder. 


2. Low-speed flow past a circular cylinder 
It will be sufficient to take both the velocity U at infinity, and the 
radius of the circle, to be unity. Then 


w=2+l1/z, ¢=1—Il1/24, z= (1—2f)-t. (14) 


The circle in the z-plane corresponds to a circle of radius unity, with 
centre at (1,0), in the ¢-plane; the z-plane outside the circle corresponds 
to the inside of the circle in the ¢-plane taken twice. In fact the flow is 
symmetrical about both the axes of x and y in the z-plane, and the values 
ofg and @ are the same in the first and third quadrants, and in the second 
and fourth quadrants. If we suppose the z-plane cut along AB and DE 
Fig. 1), the upper half of the plane corresponds to the inside of the circle 
inthe ¢-plane with a cut from the origin to (1,0) (Fig. 2a), and similarly 
ior the lower half (Fig. 26). There is a branch point at £ = 1, correspond- 
ing to the velocity at infinity, and the representation of the flow plane 
n the hodograph plane is a Riemann surface of two sheets. 

When we express z and w as series of powers of £, two series for (1—{)}, 
ne of which is minus the other, have to be taken into account even 
when |f) = g <1; near D, for example, we must take the series which 
reduces to +1 when £ = 0, and near B the series which reduces to —1. 
There are similarly two other expansions, one of which is minus the other, 
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to be taken into account when q > 1; and there are therefore four series 
representations in all. In the regions (1), (2), (3), (4) marked in Figs. 1+ 
and 2 these series are (with x! written for ['(a+1)) 





= oe = | \ 

Ww, == (n _ 3)! en 
—2)in! 
n=0 ( 2) 
oe x 

ss ; (n+ 4)(n—$)! yn r (15) 
- (—piat * ’ 

n=0 eis 
Ws = —W\ W, = —Ws, J 























% 
, a he J (3) \A 
2 2) Je rn 4) }c 
(3) , y 
Fic. 2a. Fic. 26. 


In this case it is easy to see that w,, ws, w, are the correct continua- 
tions of w, but this will now be verified by a method which is applicable 
to high-speed flow. The method is exactly analogous to the use of an 
integral of Barnes’s type to find the continuation of the hypergeometric 
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series (Whittaker and Watson, 7); in fact, apart from the factor n—1, 
w, is a hypergeometric series. 


p . — | = ! 
Consider : i 3. Naa Lt) (—2)" dy, (16) 
27rd . (—}4)! 


where the path of integration is chosen so that the poles of (v—3)!, which 
we at v = 4—n (n = 0,1, 2,...), are on the left of the path, and the poles 
of (v—1)(—v—1)!, which are at v = 0, 2, 3,..., are on the right of the path. 
We ensure convergence by taking |arg(—{)| < 7; hence we take 


c ge” 6) (0 < @ - 277) | (17) 
get) (ln <0 <0) J 

It may now be proved, exactly as in the proof of Barnes’s integral for 
the hypergeometric function, that the integral round a semicircle of 
radius V-+-4 on the right of the imaginary axis with centre at the origin, 
where V is an integer, > 0 as N +00 if g <1. The expression in (16) is 
therefore equal to minus the limit of the sum of the residues of the 
integrand on the right of the contour in (16), which is the series w, in (15). 
Qn the other hand, when g > 1, the integral round a semicircle of radius 
Y on the left of the imaginary axis with centre at the origin > 0 as N > 0, 
ind the expression in (16) is equal to the limit of the sum of the residues 
f the integrand on the left of the contour in (16), which is the series w, 
in (15) if 8 is negative and the series w, (= —w,) if @ is positive. 

Since @ is negative on the dividing lines between the regions | and 2, 
nd 3 and 4, and positive on the dividing lines between the regions 2 and 3, 
nd 4and 1, it follows that w, is the correct continuation of w,, ws (= —w) 
the correct continuation of w,, w, (= —w,) the correct continuation of ws, 


and w, the correct continuation of wy. 


3. The corresponding high-speed flow 
We take the stream-function to be the imaginary part of a function W, 
ind we start from the series corresponding to the series for w, in (15), 
namely, : 2)! 
; : + (nN \(m—#): i , 
m=)» Salt Wala dein (18) 


ra ( 5). 1: 
n=0 2 


7, 1s the value of 7 corresponding to the undisturbed velocity at infinity, 
vrhich we take to be subsonic, so that 
9; 9 ¢ 
T <7; = 95/96 = (2B+1)>. (19) 


We shall use the following properties of the functions y%,(7) (Light- 
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I. For 0 <7 <1, %,(r7) is an analytic function of n except at n = —2 
—3, —4...., where it has simple poles. The residue at n = —m jg 


—mC,, w,,(7), where 
! ! 9—2 
(a,,—1)!(m—b,,)! e—2om 


(m!)?(a,,—m—1)!(—b,,)! 2am 





vay] om 


as m > oo, a,,, b,, are defined in (10), andt 


m 


o = —7,'tanh-'7}—} log 48 
= —r;}log(rz3+1)+ 4H(rz#+ Ilog 2+ 4(7z#— Ilog B. (21) 


Il. If 5 and e are any two positive numbers, 


_ @\8841742 
V(r) = i" 7) r (22) 
_ |—t/r, 
T.—T\4 1—r/7.\} 
and $= o-+r;4tanb-{ s "J —tanh-(* ~ aie (23) 
l—r l—r 
a ] 
then as |n| > 00, ws, (7) = Vi(r)e “[1+40())], (24) 
n 


uniformly for 0 <7 <7,—e and for n in the whole complex plane if 
jn+m| > 4 for all integers m. 

tT ~ e* as t> 0 ahd s > —oo, s increases when 7 increases, and s = o 
when rt = 7,. 

Now we have already seen that in order that (18) should reduce to the 
corresponding series for low-speed flow when 7 and 7, > 0 we should like 
Si(t1) ~ Tr" as 7,>0, ie. f,(7,) ~ e-™ as 8, > —oo, where 3, is the 
value of s corresponding to t = 7,. Next we should like the behaviour 
of the series (18) near 7 = 7, to be similar to that of the first series of (15) 
near q = 1. If, for large n, f,,(7,) behaves like a function of 7, multiplied 
by e-”*1, it follows from the asymptotic formula for ¢,,(7,) that the series 
(18) will behave like the first of (15) with exp(s—s,—i@) in place of ¢, and, 
in particular, its behaviour near 7 = 7,, s = 8,, will be the same as that 
of the first of (15) near g= 1. The simplest function satisfying these 
requirements on f,(7,) is e~"*: (Lighthill, 4); so we start with 


x 


, (n—1)(n—32)! Sicaal 9x) 
W, = = r)e—M(si +00) (29) 
. pa (—4)!n! ¥n(7) 


n=0 


To find the continuation of this series. consider 


F te Si —Sitt—»— 2)! ‘i 
: ee | heme (—1)"b,(r)e-" i‘) dy, (26) 
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where 1 is e’” for positive 6 (0 <@ <2) and e-‘” for negative 
§(—27 < @ <0), and the path of integration passes to the right of the 


)!&,(r) and to the left of the poles of (»—1)(—v—1)! 
Use of the asymptotic expression for %,(7), together with a procedure 


poles of (v 


exactly similar to that for low-speed flow, shows that the integral is 
convergent, and that for 7 < 7,, in which case s < s,, the integral round 


} on the right of the imaginary axis > 0 as 


1 semicircle of radius NV 
V +00, so that (26) is equal to minus the limit of the sum of the residues 
of the integrand on the right of the path in (26), which is the series W, 
in (25). On the other hand, when 7, < + <7,, the integral round a semi- 
circle of radius N+} on the left of the imaginary axis with centre at the 
origin > 0 as N + oo, and the expression (26) is now equal to the limit of 
the sum of the residues of the integrand on the left of the path in (26). 
These residues arise partly from the poles of (v—#%)! and partly from the 
poles of ys,(7). The former give rise to the series 


, - Xo (n+4})(n—$)! . - 
W, ; S (n , ae iby (re -4Xs1 +70) (27) 
n 0 2/° ‘ 


if @ is negative, and to a series —W, if 6 is positive. The latter give rise 
toa series — L where 


— (—])™+l(m-+1)!(—m—3)! a 
—_, 8 On thm reer 


int, (—})! 
. 9m+l {C7 
* - (77 1)! m us (r)em(si+t) (28) 
2 - °: m 
— L.3.5....(2m +1) 
m 2 


Thus the continuation of the series W, is WW—Z when 6 is negative 
and —W,— L when @ is positive; the continuations into the regions corre- 
sponding to (2), (3), and (4) are therefore W,— L, —W,—2L, and —W,—L. 
It should be specially noted that the continuation of —W,—L from the 
region corresponding to (4) into the region corresponding to (1) is W,, so 
that W returns to its original value after a double circuit of the branch 
point, and ys is single-valued in the physical plane, as it is required to be. 

The resulting solution is no longer symmetrical about the axis of y, but 

¢ obtain a symmetrical solution if we start from the series Wj+L. The 
series which hold in the regions corresponding to (1), (2), (3), and (4) are 
then +L, W,, —W,—L, and —W, respectively. 


The following inequa 





ities may be deduced from results given by 
uighthill (5). When 0 <7<7,, and |n+m! > 65 for all non-negative 


ntegers m 


(r)| < Aln|*\ens|, (29) 
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where the constant A depends on 6 only. When 7, <7 < 1—e, and 
n+m| > 6 for all non-negative integers m, and n is real, 
Iyp,,(7)| < B\n|tene, (30) 


where B depends on 6 and « only; 6 and « are any two positive numbers, 
In each case the factor |n|* may be omitted if 7 is restricted to 
0 <7<7,—e or 7,+€ <7 < 1l—e, respectively, and A depends on § 
and e. 

From these inequalities it follows that the series (27) for W, converges 
not only for 7, <7 <7,, but throughout the whole of the supersonic 
region in the flow, in fact for 7,< 7 < 1—e; the series (28) for L is 
convergent everywhere, at subsonic, sonic, and supersonic speeds. We 
therefore have the continuations of W, (or Wj4+L) throughout the whole 
flow field, including the supersonic region. This is a particular case of 
a general result (Lighthill, 4). 

It is now clear that exactly the same process can be carried through 
with other values for f,(7,) which have the required behaviour as 7, > 0 
and as n > 00, and which are analytic functions of » apart from poles. 
In general, extra terms will arise from the residues at the poles of f,; if f, 
has poles at the poles of (v—1)(v—#)!(—v—1)!%,(7), double poles will 
result and the series will be more complicated. 

The same method may be used to find the continuation of more 
general series than (25), and to find the high-speed flow corresponding to 
low-speed flows other than that past a circular cylinder. Some more 
general flows have been constructed in this way by Craggs. 


4. Connexion with Lighthill’s general solution 

We may now show how the results of the preceding section are con- 
nected with Lighthill’s general solution for a subsonic region, and his 
method of obtaining series when series such as (15) are known for w. He 
constructed a solution for a general value of f,,(7,), but for our purposes 
it will be sufficient to restrict ourselves to the simple case f,,(7,) = e-"": 
the solution may then be written 


@ 


W = w()t+ > Crpn(r)em*L,, (31) 
m=2 
A 
where  gmbe gy 4 pmsl? gy (32) 
m . dt : Yi 
0 0 
and A = exp(s—s,—i6). (33) 


An arbitrary lower limit, depending on m if necessary, may be taken 
instead of zero in the integral, but this is not needed for our purposes. 
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In the first place, the connexion of this result with the method of §3 
for the particular case there considered is found by inserting into (26) the 
expansion of #,(7) in rational fractions 


x 


m & ems (r) 
/, (r vs} Jtyp m~- Tm 34 
w, (7) é | v Pa ae (34) 


(Lighthill, 5), which is the result of applying Mittag-Leffler’s theorem 
(Whittaker and Watson, 8) to e-’%/,(r). The order of summation and 
integration is then interchanged, and use made of (15) to express J. 
The manipulations may in fact be justified, but it is preferable to verify 
independently that the imaginary part of (31) is a solution of the equation 
for w. 

When A lies in the region (1) of the ¢-plane, the series (18) is obtained 
from (31) by substituting the first of (15) into J, to give 

] S (n—1)(m—$)! ndnem 


‘sp Ly —1)!n! n+m’ 
7 0 = , 


(35) 


rearrangement of the resulting double series is easily justified, and 
(34) is then applied to prove that (31) is equal to the series (25). But 
when A is in the region (2), the same procedure cannot be followed by 
substituting the second of (15) into the expression for J,,, since the path 
of integration, from 0 to A, does not lie wholly in the domain of convergence 
of the second of (15). Nevertheless it is correct that 


x « 
SS (nt4)(n—¥! (J—n)amet—n ” 
& U > 1 ’ ' = 1 +hin (36) 
ine. (—+$)in! m+s—n 


where h,, is a constant, since d/,,/dA has the correct value, A"(dw,/dl);_); 


also h,, is the constant term in the expansion of J,, in descending powers 
of \ when A is in region (2). (Similar results hold for regions (3) and (4).) 
If we again rearrange the double series and use (34), we find that for A in 
the region (2), (31) is equal to 


i 


Wt > bon Cn thn (rremer®, (37) 


m m 
where W, is the series (27). 

There are various ways of setting out the determination of the constants 
h,,; We choose that which will also enable us to arrive at the expression 
47) below. 

We find J,, in closed form in terms of w and its integrals by continued 
integration by parts. Write 

wr, (C) | w(l) dZ, wi, (Cf) = Wi» —(C) al. (38) 
i 1 
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In region (1), 


(—2)"(1—£)"-4(2m+2—L) 


w = (1—Z)'+(1—Z)-+, wi,,(f) = 3 3 (m+) ; (39) 
in region (2) | 
w= —ilt( 1) 72-4(1—{-1)-3, 
scl Qmirm-h/ -{- 1)m— *( 2m+ 2—f) 
WI» (S) = 1.3.5....(2m-+-1) (40) 


(where (1—Z)! is to be interpreted as the series that reduces to +1 when 
¢ = 0); in region (3) the values are minus those in region (1), and in 
region (4) minus those in region (2). Then by continued integration by 
parts 


I - { mora C = A™w(A)— mA" WI, (A) +... + 


m 


FY 
+(—1)"-'m(m—1)...3. 2Awi,,, (A) +(—1)™m! [ wi, (A)— wi, (0) ].. (41) 


Now w(A), wi,(A),..., wi,,(A), with A in region (2), may be expanded in a series 
of descending powers of A, involving powers which are half odd integers 
and with no constant terms; the resulting series when these terms are 
added is the series on the right in (36), and the constant term is 


ee 1 )! 


9 
— Bi (42) 
san . (2m 1) 


h, = (—1)"4'm! wi,,(0) = — 
With this value of h,,, the second term in (37) is —L, where L is given by 
(28), and the continuation of W, is W—L as before. 

We may proceed similarly in region (3) by taking the path of integration 
for I, to be a contour from the origin D in the ¢-plane (Fig. 2a), which 
encircles the point (1,0) once in the negative sense, returns to the origin 


B, and thence proceeds to A entirely in region (3). Hence with the first 
integral entirely in region (3), 


A d a=) j 
dw dw 
i | cm. df+ ym _ d 
m= | oar ir Tae 
0 0 
nn a 3) ! n+m 
- += — n—#)! = cel H,, (43) 
— (—4)!n! n+m 
=) ] 
where i = =¢ gm z at = 2, (44) 


0 
and h,, is given by (42); this last result is obtained immediately by 


continued partial integration. 
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The result is that the continuation of W, into region (3) is —W—2L, 
as before. 

When we pass from region (3) to region (4) (Fig. 26), we increase the 
value of the constant from its value in (44) by minus the value (42) found 
in passing from region (1) to region (2), so the constant becomes simply 
},,, and the continuation into region (4) is —W,—L, as before. 

The two methods have now been shown to lead to identical results. 
If we substitute from (41) into (31), invert the order of summation, and 
use the relation 


V(r) = 1+ > C,e™p,(7), (45) 


obtained by letting v > 00 in (34), we find that (with C, = 0) (31) may 
be expressed in the form 


V(r)w(A) 1 > h Cc; wb (rye m $1 +i6) 1 


x 


2 ON) 
S (—1)" : we >” m(m—1)...(m—n-+-1)C,, ,,(7)e™. (46) 
n=1 : m=n 


The second term is — L, and we are therefore led to the following form 
for W,+ L and its continuations for subsonic r: 


x . a 
: wt, (A) = ' 
V(r)w(A) + S (—1)" \n > m(m—1)...(m—n+1)C,, o,,(r)e”™. 
— é 
n=1 m=n - 
(47) 
The result may be useful for computation when A is near 1: the inner sums 
may be computed by considering the asymptotic expansion of e~"*¢, (7) 
for large n. The double series converges absolutely when 


1—A exp(2o—s—s,)—exp(s—S,). 


5. The position coordinates 

In addition to determining y (as the imaginary part of W) we have to 
letermine the position coordinates X, Y in the plane of the high-speed 
low. It is a simple matter to write down formal series for X and Y corre- 
sponding to the series W,, Wj, and L for W; %,(7) does not occur, and it 
is easy to prove that, corresponding to % = y,(7)sin(v8+e), where « is 


arbitrary and vy + 1, apart from additive constants, 
fe (a -30-1- , 
Z=X-+i¥ . jexplice EO 1) — alti — 
Vqt (] T) v+1 \2 J 
aa ane FN dato). (48) 


the formal series for X and Y are obtained by taking real and imaginary 
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parts, inserting correct values for « (« = 0 for the terms of W, and L and 
e = 47 for the terms of W,), multiplying by the multipliers and summing, 

The velocity potential may also be found if desired, since, correspond- 
ing to / = ¢,(7)sin(v8+e) (v ~ 1), 


d — — a yas(r}e0s(o0-+e). (49) 
6. Generalization of the method 
To generalize the method to apply to any initial series (13), we require 
a function F(v) with the following properties: 
(i) F(v) is analytic everywhere except at poles, which are all to the left 
of some line R(v) = constant. 
(ii) For positive integral n, F(n) = d,. 
(iii) When |v| +00 with R(v) bounded, F(v) = O(e->”') for any 8 > 0. 
(iv) When R(v) > +00, F(v) = O(e4"™) for some A. 
(v) When R(v) > —oo, F(v) = O(e-?*™) for some B. 
Then if f,(7) satisfies the conditions (i), (ii), and (iii) of ref. (4), § 2, and 
also condition (i) above, the integral 
1 +» _ , 
ami | [arene POMotedfran— Ure“ d - 
(where as in (26) —1 is e’” for @ positive and e-‘7 for @ negative) taken 
round a contour with the positive integers to the right and all other poles 
of the integrand to the left, converges; also it is equal to (13) when 
8 <8,—A and to the sum of the residues of the integrand at its poles 
not positive integers when s > s,+ B. Thus it affords a continuation of 
(13) into the region s > s,+-B, and further continuation will often be 
possible in a similar manner. 
A convenient form for F(v) is often 


(0+) 
“ 


- | o-*-tw(C) dé, (51) 
270 . 
x 
within the range of convergence of the integral, outside which its analytical 
continuation must be employed. (It is assumed that no singularity of 
w(C) lies within the contour.) (51) can be rewritten 





(0+) (20+) 
; £<e ™ LL oF ‘a 
-— w'({)df = —. C-"+1 dz, (52) 
2m Jv 27iv . 
—o Zo 


in some range of v; this reduces it to an integral in the physical plane of 
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L and the low-speed flow; unfortunately the point z,, where € = 0, must in 
nming general be inside the body. Thus it appears that the present method can 
‘spond- be applied only if the analytical properties of the low-speed potential 


inside the body are known. This is a problem intimately connected with 
that of finding a ‘generalized Laurent series’ for w(f), discussed in ref. 


(9) b (4), $4. 
It seems that no theory as general as that of ref. (4), §2, could be based 
on the present method. 
require 
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SUMMARY 


In this paper attention is confined to the steady irrotational isentropic flow of a 
perfect gas at supersonic speeds, in particular, flow in which one characteristic is 
straight. In two dimensions any characteristic along which the velocity is constant 
in magnitude and direction is straight and, conversely, if a characteristic is straight 
the velocity must be constant along it ; in three-dimensional flow with axial symmetry 
the velocity along a straight characteristic is not generally constant but varies 
according to a definite relation with the radial distance. We establish this relation 
and describe the types of flow corresponding to it. 

Notation 
r distance from the axis of symmetry 
6 angle between the direction of flow and the axis 
w magnitude of the velocity 
W non-dimensional form of the velocity = w/wyax 
a local velocity of sound 
u Mach angle 
M local Mach number 


The equation for the velocity along a straight characteristic 


D2 
BERNOULLI’S equation ha +w* = w%. (1) 
pi 
») 
may be written, replacing a by wsiny and with k? = ——.,, in the form 
on 
/ 
W = Wmax (1+F*sin®p)-!. (2) 


One condition satisfied along a characteristic making angle e with the axis 


of symmetry is : * 
' . 6 = e+up. (3) 


The choice of sign in equation (3) depends on the family to which the 
characteristic belongs. A characteristic to which the upper sign applies is 
said to be of the ‘-+-’ family, one to which the lower sign applies is of the 
‘_’ family. The same distinction is used in all following equations. 
Combining (2) and (3) we obtain 
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1+ k?sin?(e—@)}-?. (4) 
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The second condition (1) satisfied along a characteristic is 


a? dé a? dr sin @ a 
w dw+ — _ : = 0, (5) 
Sim pe COS pL r cos w sme 
which becomes, after eliminating w and @ by use of (2) and (3), 
1—k?+ 2k?sin?u 
f du 


sinysin(e+p)dr , 
1+ k*sin?u sine a 


ide = 0. (6) 
If the characteristic is straight « is constant along it. Equation (4) then 
shows that the hodograph of a straight characteristic is an ellipse circum- 
scribing the sonic circle and inscribed in the circle of maximum velocity, with 
major axis parallel to the direction of the characteristic. 

With e constant, equation (6) reduces to 


dr sin e(1 —k?+- 2k*sin2) d T{(1—k?)+(1+k?)t?} 


r sin » sin(e+){1+k*sin*y} _ “(7'+t){1+(1+k?)e} * 
if we write ¢ = tanp, T = tane. 
Hence, integrating, 
R = B-P (D+ pA DTAYN+A+k IT} yy 
[{ (is ke2)g2\Ti2exitan— (+k Milk NUTR +k ML HA ke (7) 


where R is the non-dimensional form of the distance from the axis incor- 
porating the constant of integration. Equations (2), (3), and (7) give 
respectively the magnitude and direction of the velocity and the distance 
from the axis in terms of the parameter p = tan-4t. 


Properties of flow with straight characteristics 

We note firstly that straight characteristics can exist along which the 
velocity is constant in magnitude and direction. If » and @ are constant 
equation (3) shows that « is constant. Equation (6) reduces to 


sin w sin(e+p) dr 0 
r sine 

Now dr/sine is the element of length along a characteristic and does not 

vanish and hence, if we exclude the vacuum flow corresponding to 


sin pu 0. ; , 
; sinfe+p)= 0, Le. O= 0. 


A characteristic with constant velocity must be straight and the flow at 
every point of it directed parallel to the axis of symmetry. This special 
type of straight characteristic is not represented in equation (7) because 
equation (6) has been integrated assuming sin(e+jp) = 0. 

In the more general case described by (7), the condition that R should 
be real and positive restricts the variation of t. 

On a ‘+’ characteristic ¢ > —7'; « and T must be negative if the axis 
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of symmetry is included but may have either sign otherwise. On a ‘—’ 


characteristic t < 7; « and 7 must always be positive. 

If the characteristic is to reach the axis of symmetry, R must vanish 
when ¢ takes its axial value +7. Thus the exponent of (7'+-t) in (7) 
must be positive, i.e. 








2 : 
|’ < |e) =tT, say (r= Ne for y = 1-4). (8) 
We can therefore distinguish between five general types of flow 
(a) Flow with a ‘+-’ characteristic, including the axis 
RL 
0-14 
——_—__—_ S 
- M=1 
0-1 


Fic. 1. Case (a). Flow on first ‘sheet’.+ 








ake 
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0-1F 
Le ——— 
0 M=1 
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Fic. 2. Case (a). Flow on second ‘sheet’. 
Here 0 > T > —7, t > —T; when —T <t <7, dM/dR is negative 


and d6/dR positive; when t > 7, dM/dR is positive and dé/dR negative. 
At t = «©, when sonic speed is reached, R takes the limiting value 


R = f( l + 2) Tie7il2v eT A+ KANE H1+k*)T?] 


Then, as ¢ increases from its axial value, R increases from zero to a maxi- 
mum at ¢ = 7 and then falls to the non-zero value R. 

The value ¢t = 7 is of special interest; it defines the point at which the 
characteristic reaches a limit line, beyond which potential flow cannot be 


+ Fies. 1-8. Behaviour of a vector, having magnitude of Mach number and direc- 
tion of flow, along straight characteristics. A, S, and L are the points of intersection 
of a straight characteristic with the axis of symmetry, sonic line, and limit line 
respectively. 
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continued. For if do denotes an element of the streamline at this point 
the acceleration along the streamline is given by 

dw dw dr dw . 
v w— — = w—sin 8. 

do dr do dr 
At t = 7, dr/dt = 0, while dw/dt and sin @ have non-zero values: therefore 
the acceleration is infinite at this point. 

Evidently this type of flow extends over two ‘sheets’. On the first the 
characteristic starts on a sonic line away from the axis and ends at the limit 
line, the velocity increasing until the limit line is reached; the characteristic 


t 


then returns on the second ‘sheet’ with increasing velocity to reach the 
axis of symmetry. 


(b) Flow with a ‘+-° characteristic, excluding the axis, T positive 


290¢ 











Fic. 3. Case (b). Flow on first ‘sheet’. 
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Fic. 4. Case (b). Flow on second ‘sheet’. 
Here ¢ can take all positive values; when 0 <t <7,dM/dR is positive 
and dé/dR negative; when t > 7, dM/dR is negative and d@/dR positive. 
As in case (a) the characteristic crosses from one ‘sheet’ of the flow to 
a second at a limit line. It starts on a sonic line on the first ‘sheet’. runs 
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inwards with M increasing and 6 decreasing up to the limit line. On the 
second ‘sheet’ the characteristic has no outer limit, M increasing and 
decreasing outwards from the limit line. Part of this flow is in the negative 
direction, 6 being initially > 90° at the sonic line. The change to the 
positive direction takes place on the first ‘sheet’ if « < cot-'+, on the 
second ‘sheet’ if « > cot—+. 








(c) Flow with a ‘+-’ characteristic, excluding the atis, T negative 

30 
RL 
15F 
4 premipnavedt 
M=1 
S 
0 


Fic. 5. Case (c). 
Here ¢ always exceeds 7, dM/dR is positive, and d0/dR negative. The 
characteristic starts at R = R with sonic velocity and extends indefinitely 
outwards with increasing velocity. 


(d) Flow with a ‘—’ characteristic, including the axis 








Fic. 6. Case (d). 

Here 0 <¢ < T <7,dM/dR and dé/dR are always positive. This flow 
covers the whole field with velocity increasing outwards from the axis. 
Singularities are absent. 

(e) Flow with a ‘—’ characteristic, excluding the axis 

Here 0<t <7, T >7; when t <7, dM/dR and d6/dR are positive: 
when ¢t > 7, dM/dr and d6/dR are negative. 

The flow again covers two ‘sheets’ on each of which the characteristic 
has no outer limit but is terminated on a limit line before reaching the 
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on tis TABLE I TABLE II 
and @ Case (a) T _ A; € —18-4 Case (b) T — 4; e = 18-4° 
cative Vl r) R M 6 
to the 2 I 7 281°3 1°00 108°4 
mn the *I4 T°05 260 1°03 94°7 
1°20 7 240 Ils 79°0 
74 1°58 20°8 224°0 1°58 57°7 
1°94 2 240 2°17 46°0 
1 2°15 7 260 2°50 42°I 
2°41 2 250 2°75 39°8 
2°64 ) 300 2°97 381 
2 2:89 340 3°33 360 
3°16 370 3°55 348 
TABLE ITI TABLE IV 
Case (c) T —l;« —45 Case(d) T=}; €«= 18-4 
Gf] R M 6 
I° ° 3716 0°00 
j I ) 5 or! 3°60 2°41 
1°136 8 o°2 3°76 3°02 
‘218 o"4 3°97 3°80 
I 1°275 06 410 4°27 
1e I2 1°310 od "22 63 
rl 3 4 4°93 
initely ; 1°39 1°0 4°30 4°95 
: 1*388 I°2 4°38 5°21 
I*4c 1°O 4°52 5°66 
I°407 7 2°0 4°64 6°07 
TABLE V 
> ne 
Case (e) 7 10; « 84-3 
VW 6 
' 1°03 9°6 
1°04 12°5 
1°09 Ig°O 
I°22 29°09 
=, | 1°60 45°0 
2°21 57°7 
2 2°70 63°1 
, 3°08 65°2 
is flow = 
2°32 66°7 
ie axXIs { 3°69 69°0 
axis. On the first ‘sheet’ the velocity increases as we follow the charac- 
ositive teristic inwards to the limit line; we then return on the second ‘sheet’ 
with outwardly increasing velocity. 
teristic rhe fields of flow corresponding to (a), (b), and (c) appear to be the first 
ino the with axial symmetry in which limit lines have been encountered; it is also 











364 VELOCITY ALONG A STRAIGHT CHARACTERISTIC 


interesting to note that in all three cases the limit line is always reached 
at the Mach number M = 2/,/(3—y) (= 1-581 for y = 1-4), whatever the 
slope of the characteristic. 


M=1 M=1 
L PL 
8 8 
Fic. 7. Case (e). Fic. 8. Case (e). 
Flow on first ‘sheet’. Flow on second ‘sheet’. 














Tables I-V show the variation along straight characteristics of Mach 
number and direction of flow with distance from the axis. y is taken as 1-4 


and values of 7’ are chosen so as to represent all five types of flow, namely, 
T = —}4, 4, —1, 4, 10 for (a), (b), (c), (d), (e) respectively. Figs. 1-8 


show, for each type of flow, the form of characteristic and the behaviour 
of a vector representing the magnitude and direction of the velocity 


along it. 
The author is indebted to Professor Goldstein and Dr. Meyer for theit 
helpful criticism of this paper. 
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